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