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