001:       SUBROUTINE ZHETRI( UPLO, N, A, LDA, IPIV, WORK, 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:       CHARACTER          UPLO
009:       INTEGER            INFO, LDA, N
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IPIV( * )
013:       COMPLEX*16         A( LDA, * ), WORK( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZHETRI computes the inverse of a complex Hermitian indefinite matrix
020: *  A using the factorization A = U*D*U**H or A = L*D*L**H computed by
021: *  ZHETRF.
022: *
023: *  Arguments
024: *  =========
025: *
026: *  UPLO    (input) CHARACTER*1
027: *          Specifies whether the details of the factorization are stored
028: *          as an upper or lower triangular matrix.
029: *          = 'U':  Upper triangular, form is A = U*D*U**H;
030: *          = 'L':  Lower triangular, form is A = L*D*L**H.
031: *
032: *  N       (input) INTEGER
033: *          The order of the matrix A.  N >= 0.
034: *
035: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
036: *          On entry, the block diagonal matrix D and the multipliers
037: *          used to obtain the factor U or L as computed by ZHETRF.
038: *
039: *          On exit, if INFO = 0, the (Hermitian) inverse of the original
040: *          matrix.  If UPLO = 'U', the upper triangular part of the
041: *          inverse is formed and the part of A below the diagonal is not
042: *          referenced; if UPLO = 'L' the lower triangular part of the
043: *          inverse is formed and the part of A above the diagonal is
044: *          not referenced.
045: *
046: *  LDA     (input) INTEGER
047: *          The leading dimension of the array A.  LDA >= max(1,N).
048: *
049: *  IPIV    (input) INTEGER array, dimension (N)
050: *          Details of the interchanges and the block structure of D
051: *          as determined by ZHETRF.
052: *
053: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
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, D(i,i) = 0; the matrix is singular and its
059: *               inverse could not be computed.
060: *
061: *  =====================================================================
062: *
063: *     .. Parameters ..
064:       DOUBLE PRECISION   ONE
065:       COMPLEX*16         CONE, ZERO
066:       PARAMETER          ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
067:      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
068: *     ..
069: *     .. Local Scalars ..
070:       LOGICAL            UPPER
071:       INTEGER            J, K, KP, KSTEP
072:       DOUBLE PRECISION   AK, AKP1, D, T
073:       COMPLEX*16         AKKP1, TEMP
074: *     ..
075: *     .. External Functions ..
076:       LOGICAL            LSAME
077:       COMPLEX*16         ZDOTC
078:       EXTERNAL           LSAME, ZDOTC
079: *     ..
080: *     .. External Subroutines ..
081:       EXTERNAL           XERBLA, ZCOPY, ZHEMV, ZSWAP
082: *     ..
083: *     .. Intrinsic Functions ..
084:       INTRINSIC          ABS, DBLE, DCONJG, MAX
085: *     ..
086: *     .. Executable Statements ..
087: *
088: *     Test the input parameters.
089: *
090:       INFO = 0
091:       UPPER = LSAME( UPLO, 'U' )
092:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
093:          INFO = -1
094:       ELSE IF( N.LT.0 ) THEN
095:          INFO = -2
096:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
097:          INFO = -4
098:       END IF
099:       IF( INFO.NE.0 ) THEN
100:          CALL XERBLA( 'ZHETRI', -INFO )
101:          RETURN
102:       END IF
103: *
104: *     Quick return if possible
105: *
106:       IF( N.EQ.0 )
107:      $   RETURN
108: *
109: *     Check that the diagonal matrix D is nonsingular.
110: *
111:       IF( UPPER ) THEN
112: *
113: *        Upper triangular storage: examine D from bottom to top
114: *
115:          DO 10 INFO = N, 1, -1
116:             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
117:      $         RETURN
118:    10    CONTINUE
119:       ELSE
120: *
121: *        Lower triangular storage: examine D from top to bottom.
122: *
123:          DO 20 INFO = 1, N
124:             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
125:      $         RETURN
126:    20    CONTINUE
127:       END IF
128:       INFO = 0
129: *
130:       IF( UPPER ) THEN
131: *
132: *        Compute inv(A) from the factorization A = U*D*U'.
133: *
134: *        K is the main loop index, increasing from 1 to N in steps of
135: *        1 or 2, depending on the size of the diagonal blocks.
136: *
137:          K = 1
138:    30    CONTINUE
139: *
140: *        If K > N, exit from loop.
141: *
142:          IF( K.GT.N )
143:      $      GO TO 50
144: *
145:          IF( IPIV( K ).GT.0 ) THEN
146: *
147: *           1 x 1 diagonal block
148: *
149: *           Invert the diagonal block.
150: *
151:             A( K, K ) = ONE / DBLE( A( K, K ) )
152: *
153: *           Compute column K of the inverse.
154: *
155:             IF( K.GT.1 ) THEN
156:                CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
157:                CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
158:      $                     A( 1, K ), 1 )
159:                A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
160:      $                     K ), 1 ) )
161:             END IF
162:             KSTEP = 1
163:          ELSE
164: *
165: *           2 x 2 diagonal block
166: *
167: *           Invert the diagonal block.
168: *
169:             T = ABS( A( K, K+1 ) )
170:             AK = DBLE( A( K, K ) ) / T
171:             AKP1 = DBLE( A( K+1, K+1 ) ) / T
172:             AKKP1 = A( K, K+1 ) / T
173:             D = T*( AK*AKP1-ONE )
174:             A( K, K ) = AKP1 / D
175:             A( K+1, K+1 ) = AK / D
176:             A( K, K+1 ) = -AKKP1 / D
177: *
178: *           Compute columns K and K+1 of the inverse.
179: *
180:             IF( K.GT.1 ) THEN
181:                CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
182:                CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
183:      $                     A( 1, K ), 1 )
184:                A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
185:      $                     K ), 1 ) )
186:                A( K, K+1 ) = A( K, K+1 ) -
187:      $                       ZDOTC( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
188:                CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
189:                CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
190:      $                     A( 1, K+1 ), 1 )
191:                A( K+1, K+1 ) = A( K+1, K+1 ) -
192:      $                         DBLE( ZDOTC( K-1, WORK, 1, A( 1, K+1 ),
193:      $                         1 ) )
194:             END IF
195:             KSTEP = 2
196:          END IF
197: *
198:          KP = ABS( IPIV( K ) )
199:          IF( KP.NE.K ) THEN
200: *
201: *           Interchange rows and columns K and KP in the leading
202: *           submatrix A(1:k+1,1:k+1)
203: *
204:             CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
205:             DO 40 J = KP + 1, K - 1
206:                TEMP = DCONJG( A( J, K ) )
207:                A( J, K ) = DCONJG( A( KP, J ) )
208:                A( KP, J ) = TEMP
209:    40       CONTINUE
210:             A( KP, K ) = DCONJG( A( KP, K ) )
211:             TEMP = A( K, K )
212:             A( K, K ) = A( KP, KP )
213:             A( KP, KP ) = TEMP
214:             IF( KSTEP.EQ.2 ) THEN
215:                TEMP = A( K, K+1 )
216:                A( K, K+1 ) = A( KP, K+1 )
217:                A( KP, K+1 ) = TEMP
218:             END IF
219:          END IF
220: *
221:          K = K + KSTEP
222:          GO TO 30
223:    50    CONTINUE
224: *
225:       ELSE
226: *
227: *        Compute inv(A) from the factorization A = L*D*L'.
228: *
229: *        K is the main loop index, increasing from 1 to N in steps of
230: *        1 or 2, depending on the size of the diagonal blocks.
231: *
232:          K = N
233:    60    CONTINUE
234: *
235: *        If K < 1, exit from loop.
236: *
237:          IF( K.LT.1 )
238:      $      GO TO 80
239: *
240:          IF( IPIV( K ).GT.0 ) THEN
241: *
242: *           1 x 1 diagonal block
243: *
244: *           Invert the diagonal block.
245: *
246:             A( K, K ) = ONE / DBLE( A( K, K ) )
247: *
248: *           Compute column K of the inverse.
249: *
250:             IF( K.LT.N ) THEN
251:                CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
252:                CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
253:      $                     1, ZERO, A( K+1, K ), 1 )
254:                A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
255:      $                     A( K+1, K ), 1 ) )
256:             END IF
257:             KSTEP = 1
258:          ELSE
259: *
260: *           2 x 2 diagonal block
261: *
262: *           Invert the diagonal block.
263: *
264:             T = ABS( A( K, K-1 ) )
265:             AK = DBLE( A( K-1, K-1 ) ) / T
266:             AKP1 = DBLE( A( K, K ) ) / T
267:             AKKP1 = A( K, K-1 ) / T
268:             D = T*( AK*AKP1-ONE )
269:             A( K-1, K-1 ) = AKP1 / D
270:             A( K, K ) = AK / D
271:             A( K, K-1 ) = -AKKP1 / D
272: *
273: *           Compute columns K-1 and K of the inverse.
274: *
275:             IF( K.LT.N ) THEN
276:                CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
277:                CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
278:      $                     1, ZERO, A( K+1, K ), 1 )
279:                A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
280:      $                     A( K+1, K ), 1 ) )
281:                A( K, K-1 ) = A( K, K-1 ) -
282:      $                       ZDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
283:      $                       1 )
284:                CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
285:                CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
286:      $                     1, ZERO, A( K+1, K-1 ), 1 )
287:                A( K-1, K-1 ) = A( K-1, K-1 ) -
288:      $                         DBLE( ZDOTC( N-K, WORK, 1, A( K+1, K-1 ),
289:      $                         1 ) )
290:             END IF
291:             KSTEP = 2
292:          END IF
293: *
294:          KP = ABS( IPIV( K ) )
295:          IF( KP.NE.K ) THEN
296: *
297: *           Interchange rows and columns K and KP in the trailing
298: *           submatrix A(k-1:n,k-1:n)
299: *
300:             IF( KP.LT.N )
301:      $         CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
302:             DO 70 J = K + 1, KP - 1
303:                TEMP = DCONJG( A( J, K ) )
304:                A( J, K ) = DCONJG( A( KP, J ) )
305:                A( KP, J ) = TEMP
306:    70       CONTINUE
307:             A( KP, K ) = DCONJG( A( KP, K ) )
308:             TEMP = A( K, K )
309:             A( K, K ) = A( KP, KP )
310:             A( KP, KP ) = TEMP
311:             IF( KSTEP.EQ.2 ) THEN
312:                TEMP = A( K, K-1 )
313:                A( K, K-1 ) = A( KP, K-1 )
314:                A( KP, K-1 ) = TEMP
315:             END IF
316:          END IF
317: *
318:          K = K - KSTEP
319:          GO TO 60
320:    80    CONTINUE
321:       END IF
322: *
323:       RETURN
324: *
325: *     End of ZHETRI
326: *
327:       END
328: