001:       SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            INFO, LDA, LWORK, N
009: *     ..
010: *     .. Array Arguments ..
011:       INTEGER            IPIV( * )
012:       DOUBLE PRECISION   A( LDA, * ), WORK( * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  DGETRI computes the inverse of a matrix using the LU factorization
019: *  computed by DGETRF.
020: *
021: *  This method inverts U and then computes inv(A) by solving the system
022: *  inv(A)*L = inv(U) for inv(A).
023: *
024: *  Arguments
025: *  =========
026: *
027: *  N       (input) INTEGER
028: *          The order of the matrix A.  N >= 0.
029: *
030: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
031: *          On entry, the factors L and U from the factorization
032: *          A = P*L*U as computed by DGETRF.
033: *          On exit, if INFO = 0, the inverse of the original matrix A.
034: *
035: *  LDA     (input) INTEGER
036: *          The leading dimension of the array A.  LDA >= max(1,N).
037: *
038: *  IPIV    (input) INTEGER array, dimension (N)
039: *          The pivot indices from DGETRF; for 1<=i<=N, row i of the
040: *          matrix was interchanged with row IPIV(i).
041: *
042: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
043: *          On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
044: *
045: *  LWORK   (input) INTEGER
046: *          The dimension of the array WORK.  LWORK >= max(1,N).
047: *          For optimal performance LWORK >= N*NB, where NB is
048: *          the optimal blocksize returned by ILAENV.
049: *
050: *          If LWORK = -1, then a workspace query is assumed; the routine
051: *          only calculates the optimal size of the WORK array, returns
052: *          this value as the first entry of the WORK array, and no error
053: *          message related to LWORK is issued by XERBLA.
054: *
055: *  INFO    (output) INTEGER
056: *          = 0:  successful exit
057: *          < 0:  if INFO = -i, the i-th argument had an illegal value
058: *          > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is
059: *                singular and its inverse could not be computed.
060: *
061: *  =====================================================================
062: *
063: *     .. Parameters ..
064:       DOUBLE PRECISION   ZERO, ONE
065:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
066: *     ..
067: *     .. Local Scalars ..
068:       LOGICAL            LQUERY
069:       INTEGER            I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
070:      $                   NBMIN, NN
071: *     ..
072: *     .. External Functions ..
073:       INTEGER            ILAENV
074:       EXTERNAL           ILAENV
075: *     ..
076: *     .. External Subroutines ..
077:       EXTERNAL           DGEMM, DGEMV, DSWAP, DTRSM, DTRTRI, XERBLA
078: *     ..
079: *     .. Intrinsic Functions ..
080:       INTRINSIC          MAX, MIN
081: *     ..
082: *     .. Executable Statements ..
083: *
084: *     Test the input parameters.
085: *
086:       INFO = 0
087:       NB = ILAENV( 1, 'DGETRI', ' ', N, -1, -1, -1 )
088:       LWKOPT = N*NB
089:       WORK( 1 ) = LWKOPT
090:       LQUERY = ( LWORK.EQ.-1 )
091:       IF( N.LT.0 ) THEN
092:          INFO = -1
093:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
094:          INFO = -3
095:       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
096:          INFO = -6
097:       END IF
098:       IF( INFO.NE.0 ) THEN
099:          CALL XERBLA( 'DGETRI', -INFO )
100:          RETURN
101:       ELSE IF( LQUERY ) THEN
102:          RETURN
103:       END IF
104: *
105: *     Quick return if possible
106: *
107:       IF( N.EQ.0 )
108:      $   RETURN
109: *
110: *     Form inv(U).  If INFO > 0 from DTRTRI, then U is singular,
111: *     and the inverse is not computed.
112: *
113:       CALL DTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
114:       IF( INFO.GT.0 )
115:      $   RETURN
116: *
117:       NBMIN = 2
118:       LDWORK = N
119:       IF( NB.GT.1 .AND. NB.LT.N ) THEN
120:          IWS = MAX( LDWORK*NB, 1 )
121:          IF( LWORK.LT.IWS ) THEN
122:             NB = LWORK / LDWORK
123:             NBMIN = MAX( 2, ILAENV( 2, 'DGETRI', ' ', N, -1, -1, -1 ) )
124:          END IF
125:       ELSE
126:          IWS = N
127:       END IF
128: *
129: *     Solve the equation inv(A)*L = inv(U) for inv(A).
130: *
131:       IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
132: *
133: *        Use unblocked code.
134: *
135:          DO 20 J = N, 1, -1
136: *
137: *           Copy current column of L to WORK and replace with zeros.
138: *
139:             DO 10 I = J + 1, N
140:                WORK( I ) = A( I, J )
141:                A( I, J ) = ZERO
142:    10       CONTINUE
143: *
144: *           Compute current column of inv(A).
145: *
146:             IF( J.LT.N )
147:      $         CALL DGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
148:      $                     LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
149:    20    CONTINUE
150:       ELSE
151: *
152: *        Use blocked code.
153: *
154:          NN = ( ( N-1 ) / NB )*NB + 1
155:          DO 50 J = NN, 1, -NB
156:             JB = MIN( NB, N-J+1 )
157: *
158: *           Copy current block column of L to WORK and replace with
159: *           zeros.
160: *
161:             DO 40 JJ = J, J + JB - 1
162:                DO 30 I = JJ + 1, N
163:                   WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
164:                   A( I, JJ ) = ZERO
165:    30          CONTINUE
166:    40       CONTINUE
167: *
168: *           Compute current block column of inv(A).
169: *
170:             IF( J+JB.LE.N )
171:      $         CALL DGEMM( 'No transpose', 'No transpose', N, JB,
172:      $                     N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
173:      $                     WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
174:             CALL DTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
175:      $                  ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
176:    50    CONTINUE
177:       END IF
178: *
179: *     Apply column interchanges.
180: *
181:       DO 60 J = N - 1, 1, -1
182:          JP = IPIV( J )
183:          IF( JP.NE.J )
184:      $      CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
185:    60 CONTINUE
186: *
187:       WORK( 1 ) = IWS
188:       RETURN
189: *
190: *     End of DGETRI
191: *
192:       END
193: