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