1 SUBROUTINE ZHPTRI( UPLO, N, AP, IPIV, WORK, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 AP( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZHPTRI computes the inverse of a complex Hermitian indefinite matrix
21 * A in packed storage using the factorization A = U*D*U**H or
22 * A = L*D*L**H computed by ZHPTRF.
23 *
24 * Arguments
25 * =========
26 *
27 * UPLO (input) CHARACTER*1
28 * Specifies whether the details of the factorization are stored
29 * as an upper or lower triangular matrix.
30 * = 'U': Upper triangular, form is A = U*D*U**H;
31 * = 'L': Lower triangular, form is A = L*D*L**H.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
37 * On entry, the block diagonal matrix D and the multipliers
38 * used to obtain the factor U or L as computed by ZHPTRF,
39 * stored as a packed triangular matrix.
40 *
41 * On exit, if INFO = 0, the (Hermitian) inverse of the original
42 * matrix, stored as a packed triangular matrix. The j-th column
43 * of inv(A) is stored in the array AP as follows:
44 * if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
45 * if UPLO = 'L',
46 * AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
47 *
48 * IPIV (input) INTEGER array, dimension (N)
49 * Details of the interchanges and the block structure of D
50 * as determined by ZHPTRF.
51 *
52 * WORK (workspace) COMPLEX*16 array, dimension (N)
53 *
54 * INFO (output) INTEGER
55 * = 0: successful exit
56 * < 0: if INFO = -i, the i-th argument had an illegal value
57 * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
58 * inverse could not be computed.
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ONE
64 COMPLEX*16 CONE, ZERO
65 PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
66 $ ZERO = ( 0.0D+0, 0.0D+0 ) )
67 * ..
68 * .. Local Scalars ..
69 LOGICAL UPPER
70 INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
71 DOUBLE PRECISION AK, AKP1, D, T
72 COMPLEX*16 AKKP1, TEMP
73 * ..
74 * .. External Functions ..
75 LOGICAL LSAME
76 COMPLEX*16 ZDOTC
77 EXTERNAL LSAME, ZDOTC
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL XERBLA, ZCOPY, ZHPMV, ZSWAP
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC ABS, DBLE, DCONJG
84 * ..
85 * .. Executable Statements ..
86 *
87 * Test the input parameters.
88 *
89 INFO = 0
90 UPPER = LSAME( UPLO, 'U' )
91 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
92 INFO = -1
93 ELSE IF( N.LT.0 ) THEN
94 INFO = -2
95 END IF
96 IF( INFO.NE.0 ) THEN
97 CALL XERBLA( 'ZHPTRI', -INFO )
98 RETURN
99 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**H.
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**H.
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
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 AP( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZHPTRI computes the inverse of a complex Hermitian indefinite matrix
21 * A in packed storage using the factorization A = U*D*U**H or
22 * A = L*D*L**H computed by ZHPTRF.
23 *
24 * Arguments
25 * =========
26 *
27 * UPLO (input) CHARACTER*1
28 * Specifies whether the details of the factorization are stored
29 * as an upper or lower triangular matrix.
30 * = 'U': Upper triangular, form is A = U*D*U**H;
31 * = 'L': Lower triangular, form is A = L*D*L**H.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
37 * On entry, the block diagonal matrix D and the multipliers
38 * used to obtain the factor U or L as computed by ZHPTRF,
39 * stored as a packed triangular matrix.
40 *
41 * On exit, if INFO = 0, the (Hermitian) inverse of the original
42 * matrix, stored as a packed triangular matrix. The j-th column
43 * of inv(A) is stored in the array AP as follows:
44 * if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
45 * if UPLO = 'L',
46 * AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
47 *
48 * IPIV (input) INTEGER array, dimension (N)
49 * Details of the interchanges and the block structure of D
50 * as determined by ZHPTRF.
51 *
52 * WORK (workspace) COMPLEX*16 array, dimension (N)
53 *
54 * INFO (output) INTEGER
55 * = 0: successful exit
56 * < 0: if INFO = -i, the i-th argument had an illegal value
57 * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
58 * inverse could not be computed.
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ONE
64 COMPLEX*16 CONE, ZERO
65 PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
66 $ ZERO = ( 0.0D+0, 0.0D+0 ) )
67 * ..
68 * .. Local Scalars ..
69 LOGICAL UPPER
70 INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
71 DOUBLE PRECISION AK, AKP1, D, T
72 COMPLEX*16 AKKP1, TEMP
73 * ..
74 * .. External Functions ..
75 LOGICAL LSAME
76 COMPLEX*16 ZDOTC
77 EXTERNAL LSAME, ZDOTC
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL XERBLA, ZCOPY, ZHPMV, ZSWAP
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC ABS, DBLE, DCONJG
84 * ..
85 * .. Executable Statements ..
86 *
87 * Test the input parameters.
88 *
89 INFO = 0
90 UPPER = LSAME( UPLO, 'U' )
91 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
92 INFO = -1
93 ELSE IF( N.LT.0 ) THEN
94 INFO = -2
95 END IF
96 IF( INFO.NE.0 ) THEN
97 CALL XERBLA( 'ZHPTRI', -INFO )
98 RETURN
99 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**H.
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**H.
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