1 SUBROUTINE ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, 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 * -- Written by Julie Langou of the Univ. of TN --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER UPLO
12 INTEGER INFO, LDA, N, NB
13 * ..
14 * .. Array Arguments ..
15 INTEGER IPIV( * )
16 COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix
23 * A using the factorization A = U*D*U**H or A = L*D*L**H computed by
24 * ZHETRF.
25 *
26 * Arguments
27 * =========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the details of the factorization are stored
31 * as an upper or lower triangular matrix.
32 * = 'U': Upper triangular, form is A = U*D*U**H;
33 * = 'L': Lower triangular, form is A = L*D*L**H.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
39 * On entry, the NNB diagonal matrix D and the multipliers
40 * used to obtain the factor U or L as computed by ZHETRF.
41 *
42 * On exit, if INFO = 0, the (symmetric) inverse of the original
43 * matrix. If UPLO = 'U', the upper triangular part of the
44 * inverse is formed and the part of A below the diagonal is not
45 * referenced; if UPLO = 'L' the lower triangular part of the
46 * inverse is formed and the part of A above the diagonal is
47 * not referenced.
48 *
49 * LDA (input) INTEGER
50 * The leading dimension of the array A. LDA >= max(1,N).
51 *
52 * IPIV (input) INTEGER array, dimension (N)
53 * Details of the interchanges and the NNB structure of D
54 * as determined by ZHETRF.
55 *
56 * WORK (workspace) COMPLEX*16 array, dimension (N+NNB+1,NNB+3)
57 *
58 * NB (input) INTEGER
59 * Block size
60 *
61 * INFO (output) INTEGER
62 * = 0: successful exit
63 * < 0: if INFO = -i, the i-th argument had an illegal value
64 * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
65 * inverse could not be computed.
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70 REAL ONE
71 COMPLEX*16 CONE, ZERO
72 PARAMETER ( ONE = 1.0D+0,
73 $ CONE = ( 1.0D+0, 0.0D+0 ),
74 $ ZERO = ( 0.0D+0, 0.0D+0 ) )
75 * ..
76 * .. Local Scalars ..
77 LOGICAL UPPER
78 INTEGER I, IINFO, IP, K, CUT, NNB
79 INTEGER COUNT
80 INTEGER J, U11, INVD
81
82 COMPLEX*16 AK, AKKP1, AKP1, D, T
83 COMPLEX*16 U01_I_J, U01_IP1_J
84 COMPLEX*16 U11_I_J, U11_IP1_J
85 * ..
86 * .. External Functions ..
87 LOGICAL LSAME
88 EXTERNAL LSAME
89 * ..
90 * .. External Subroutines ..
91 EXTERNAL ZSYCONV, XERBLA, ZTRTRI
92 EXTERNAL ZGEMM, ZTRMM, ZHESWAPR
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC MAX
96 * ..
97 * .. Executable Statements ..
98 *
99 * Test the input parameters.
100 *
101 INFO = 0
102 UPPER = LSAME( UPLO, 'U' )
103 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
104 INFO = -1
105 ELSE IF( N.LT.0 ) THEN
106 INFO = -2
107 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
108 INFO = -4
109 END IF
110 *
111 * Quick return if possible
112 *
113 *
114 IF( INFO.NE.0 ) THEN
115 CALL XERBLA( 'ZHETRI2X', -INFO )
116 RETURN
117 END IF
118 IF( N.EQ.0 )
119 $ RETURN
120 *
121 * Convert A
122 * Workspace got Non-diag elements of D
123 *
124 CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
125 *
126 * Check that the diagonal matrix D is nonsingular.
127 *
128 IF( UPPER ) THEN
129 *
130 * Upper triangular storage: examine D from bottom to top
131 *
132 DO INFO = N, 1, -1
133 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
134 $ RETURN
135 END DO
136 ELSE
137 *
138 * Lower triangular storage: examine D from top to bottom.
139 *
140 DO INFO = 1, N
141 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
142 $ RETURN
143 END DO
144 END IF
145 INFO = 0
146 *
147 * Splitting Workspace
148 * U01 is a block (N,NB+1)
149 * The first element of U01 is in WORK(1,1)
150 * U11 is a block (NB+1,NB+1)
151 * The first element of U11 is in WORK(N+1,1)
152 U11 = N
153 * INVD is a block (N,2)
154 * The first element of INVD is in WORK(1,INVD)
155 INVD = NB+2
156
157 IF( UPPER ) THEN
158 *
159 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
160 *
161 CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
162 *
163 * inv(D) and inv(D)*inv(U)
164 *
165 K=1
166 DO WHILE ( K .LE. N )
167 IF( IPIV( K ).GT.0 ) THEN
168 * 1 x 1 diagonal NNB
169 WORK(K,INVD) = ONE / REAL ( A( K, K ) )
170 WORK(K,INVD+1) = 0
171 K=K+1
172 ELSE
173 * 2 x 2 diagonal NNB
174 T = ABS ( WORK(K+1,1) )
175 AK = REAL ( A( K, K ) ) / T
176 AKP1 = REAL ( A( K+1, K+1 ) ) / T
177 AKKP1 = WORK(K+1,1) / T
178 D = T*( AK*AKP1-ONE )
179 WORK(K,INVD) = AKP1 / D
180 WORK(K+1,INVD+1) = AK / D
181 WORK(K,INVD+1) = -AKKP1 / D
182 WORK(K+1,INVD) = DCONJG (WORK(K,INVD+1) )
183 K=K+2
184 END IF
185 END DO
186 *
187 * inv(U**H) = (inv(U))**H
188 *
189 * inv(U**H)*inv(D)*inv(U)
190 *
191 CUT=N
192 DO WHILE (CUT .GT. 0)
193 NNB=NB
194 IF (CUT .LE. NNB) THEN
195 NNB=CUT
196 ELSE
197 COUNT = 0
198 * count negative elements,
199 DO I=CUT+1-NNB,CUT
200 IF (IPIV(I) .LT. 0) COUNT=COUNT+1
201 END DO
202 * need a even number for a clear cut
203 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
204 END IF
205
206 CUT=CUT-NNB
207 *
208 * U01 Block
209 *
210 DO I=1,CUT
211 DO J=1,NNB
212 WORK(I,J)=A(I,CUT+J)
213 END DO
214 END DO
215 *
216 * U11 Block
217 *
218 DO I=1,NNB
219 WORK(U11+I,I)=CONE
220 DO J=1,I-1
221 WORK(U11+I,J)=ZERO
222 END DO
223 DO J=I+1,NNB
224 WORK(U11+I,J)=A(CUT+I,CUT+J)
225 END DO
226 END DO
227 *
228 * invD*U01
229 *
230 I=1
231 DO WHILE (I .LE. CUT)
232 IF (IPIV(I) > 0) THEN
233 DO J=1,NNB
234 WORK(I,J)=WORK(I,INVD)*WORK(I,J)
235 END DO
236 I=I+1
237 ELSE
238 DO J=1,NNB
239 U01_I_J = WORK(I,J)
240 U01_IP1_J = WORK(I+1,J)
241 WORK(I,J)=WORK(I,INVD)*U01_I_J+
242 $ WORK(I,INVD+1)*U01_IP1_J
243 WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
244 $ WORK(I+1,INVD+1)*U01_IP1_J
245 END DO
246 I=I+2
247 END IF
248 END DO
249 *
250 * invD1*U11
251 *
252 I=1
253 DO WHILE (I .LE. NNB)
254 IF (IPIV(CUT+I) > 0) THEN
255 DO J=I,NNB
256 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
257 END DO
258 I=I+1
259 ELSE
260 DO J=I,NNB
261 U11_I_J = WORK(U11+I,J)
262 U11_IP1_J = WORK(U11+I+1,J)
263 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
264 $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
265 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
266 $ WORK(CUT+I+1,INVD+1)*U11_IP1_J
267 END DO
268 I=I+2
269 END IF
270 END DO
271 *
272 * U11**H*invD1*U11->U11
273 *
274 CALL ZTRMM('L','U','C','U',NNB, NNB,
275 $ CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
276 *
277 DO I=1,NNB
278 DO J=I,NNB
279 A(CUT+I,CUT+J)=WORK(U11+I,J)
280 END DO
281 END DO
282 *
283 * U01**H*invD*U01->A(CUT+I,CUT+J)
284 *
285 CALL ZGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA,
286 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
287 *
288 * U11 = U11**H*invD1*U11 + U01**H*invD*U01
289 *
290 DO I=1,NNB
291 DO J=I,NNB
292 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
293 END DO
294 END DO
295 *
296 * U01 = U00**H*invD0*U01
297 *
298 CALL ZTRMM('L',UPLO,'C','U',CUT, NNB,
299 $ CONE,A,LDA,WORK,N+NB+1)
300
301 *
302 * Update U01
303 *
304 DO I=1,CUT
305 DO J=1,NNB
306 A(I,CUT+J)=WORK(I,J)
307 END DO
308 END DO
309 *
310 * Next Block
311 *
312 END DO
313 *
314 * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
315 *
316 I=1
317 DO WHILE ( I .LE. N )
318 IF( IPIV(I) .GT. 0 ) THEN
319 IP=IPIV(I)
320 IF (I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP )
321 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I )
322 ELSE
323 IP=-IPIV(I)
324 I=I+1
325 IF ( (I-1) .LT. IP)
326 $ CALL ZHESWAPR( UPLO, N, A, LDA, I-1 ,IP )
327 IF ( (I-1) .GT. IP)
328 $ CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I-1 )
329 ENDIF
330 I=I+1
331 END DO
332 ELSE
333 *
334 * LOWER...
335 *
336 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
337 *
338 CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
339 *
340 * inv(D) and inv(D)*inv(U)
341 *
342 K=N
343 DO WHILE ( K .GE. 1 )
344 IF( IPIV( K ).GT.0 ) THEN
345 * 1 x 1 diagonal NNB
346 WORK(K,INVD) = ONE / REAL ( A( K, K ) )
347 WORK(K,INVD+1) = 0
348 K=K-1
349 ELSE
350 * 2 x 2 diagonal NNB
351 T = ABS ( WORK(K-1,1) )
352 AK = REAL ( A( K-1, K-1 ) ) / T
353 AKP1 = REAL ( A( K, K ) ) / T
354 AKKP1 = WORK(K-1,1) / T
355 D = T*( AK*AKP1-ONE )
356 WORK(K-1,INVD) = AKP1 / D
357 WORK(K,INVD) = AK / D
358 WORK(K,INVD+1) = -AKKP1 / D
359 WORK(K-1,INVD+1) = DCONJG (WORK(K,INVD+1) )
360 K=K-2
361 END IF
362 END DO
363 *
364 * inv(U**H) = (inv(U))**H
365 *
366 * inv(U**H)*inv(D)*inv(U)
367 *
368 CUT=0
369 DO WHILE (CUT .LT. N)
370 NNB=NB
371 IF (CUT + NNB .GE. N) THEN
372 NNB=N-CUT
373 ELSE
374 COUNT = 0
375 * count negative elements,
376 DO I=CUT+1,CUT+NNB
377 IF (IPIV(I) .LT. 0) COUNT=COUNT+1
378 END DO
379 * need a even number for a clear cut
380 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
381 END IF
382 * L21 Block
383 DO I=1,N-CUT-NNB
384 DO J=1,NNB
385 WORK(I,J)=A(CUT+NNB+I,CUT+J)
386 END DO
387 END DO
388 * L11 Block
389 DO I=1,NNB
390 WORK(U11+I,I)=CONE
391 DO J=I+1,NNB
392 WORK(U11+I,J)=ZERO
393 END DO
394 DO J=1,I-1
395 WORK(U11+I,J)=A(CUT+I,CUT+J)
396 END DO
397 END DO
398 *
399 * invD*L21
400 *
401 I=N-CUT-NNB
402 DO WHILE (I .GE. 1)
403 IF (IPIV(CUT+NNB+I) > 0) THEN
404 DO J=1,NNB
405 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
406 END DO
407 I=I-1
408 ELSE
409 DO J=1,NNB
410 U01_I_J = WORK(I,J)
411 U01_IP1_J = WORK(I-1,J)
412 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
413 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
414 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
415 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
416 END DO
417 I=I-2
418 END IF
419 END DO
420 *
421 * invD1*L11
422 *
423 I=NNB
424 DO WHILE (I .GE. 1)
425 IF (IPIV(CUT+I) > 0) THEN
426 DO J=1,NNB
427 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
428 END DO
429 I=I-1
430 ELSE
431 DO J=1,NNB
432 U11_I_J = WORK(U11+I,J)
433 U11_IP1_J = WORK(U11+I-1,J)
434 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
435 $ WORK(CUT+I,INVD+1)*U11_IP1_J
436 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
437 $ WORK(CUT+I-1,INVD)*U11_IP1_J
438 END DO
439 I=I-2
440 END IF
441 END DO
442 *
443 * L11**H*invD1*L11->L11
444 *
445 CALL ZTRMM('L',UPLO,'C','U',NNB, NNB,
446 $ CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
447 *
448 DO I=1,NNB
449 DO J=1,I
450 A(CUT+I,CUT+J)=WORK(U11+I,J)
451 END DO
452 END DO
453 *
454 IF ( (CUT+NNB) .LT. N ) THEN
455 *
456 * L21**H*invD2*L21->A(CUT+I,CUT+J)
457 *
458 CALL ZGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1)
459 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
460
461 *
462 * L11 = L11**H*invD1*L11 + U01**H*invD*U01
463 *
464 DO I=1,NNB
465 DO J=1,I
466 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
467 END DO
468 END DO
469 *
470 * L01 = L22**H*invD2*L21
471 *
472 CALL ZTRMM('L',UPLO,'C','U', N-NNB-CUT, NNB,
473 $ CONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
474
475 * Update L21
476 DO I=1,N-CUT-NNB
477 DO J=1,NNB
478 A(CUT+NNB+I,CUT+J)=WORK(I,J)
479 END DO
480 END DO
481 ELSE
482 *
483 * L11 = L11**H*invD1*L11
484 *
485 DO I=1,NNB
486 DO J=1,I
487 A(CUT+I,CUT+J)=WORK(U11+I,J)
488 END DO
489 END DO
490 END IF
491 *
492 * Next Block
493 *
494 CUT=CUT+NNB
495 END DO
496 *
497 * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
498 *
499 I=N
500 DO WHILE ( I .GE. 1 )
501 IF( IPIV(I) .GT. 0 ) THEN
502 IP=IPIV(I)
503 IF (I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP )
504 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I )
505 ELSE
506 IP=-IPIV(I)
507 IF ( I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP )
508 IF ( I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I )
509 I=I-1
510 ENDIF
511 I=I-1
512 END DO
513 END IF
514 *
515 RETURN
516 *
517 * End of ZHETRI2X
518 *
519 END
520
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 * -- Written by Julie Langou of the Univ. of TN --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER UPLO
12 INTEGER INFO, LDA, N, NB
13 * ..
14 * .. Array Arguments ..
15 INTEGER IPIV( * )
16 COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix
23 * A using the factorization A = U*D*U**H or A = L*D*L**H computed by
24 * ZHETRF.
25 *
26 * Arguments
27 * =========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the details of the factorization are stored
31 * as an upper or lower triangular matrix.
32 * = 'U': Upper triangular, form is A = U*D*U**H;
33 * = 'L': Lower triangular, form is A = L*D*L**H.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
39 * On entry, the NNB diagonal matrix D and the multipliers
40 * used to obtain the factor U or L as computed by ZHETRF.
41 *
42 * On exit, if INFO = 0, the (symmetric) inverse of the original
43 * matrix. If UPLO = 'U', the upper triangular part of the
44 * inverse is formed and the part of A below the diagonal is not
45 * referenced; if UPLO = 'L' the lower triangular part of the
46 * inverse is formed and the part of A above the diagonal is
47 * not referenced.
48 *
49 * LDA (input) INTEGER
50 * The leading dimension of the array A. LDA >= max(1,N).
51 *
52 * IPIV (input) INTEGER array, dimension (N)
53 * Details of the interchanges and the NNB structure of D
54 * as determined by ZHETRF.
55 *
56 * WORK (workspace) COMPLEX*16 array, dimension (N+NNB+1,NNB+3)
57 *
58 * NB (input) INTEGER
59 * Block size
60 *
61 * INFO (output) INTEGER
62 * = 0: successful exit
63 * < 0: if INFO = -i, the i-th argument had an illegal value
64 * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
65 * inverse could not be computed.
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70 REAL ONE
71 COMPLEX*16 CONE, ZERO
72 PARAMETER ( ONE = 1.0D+0,
73 $ CONE = ( 1.0D+0, 0.0D+0 ),
74 $ ZERO = ( 0.0D+0, 0.0D+0 ) )
75 * ..
76 * .. Local Scalars ..
77 LOGICAL UPPER
78 INTEGER I, IINFO, IP, K, CUT, NNB
79 INTEGER COUNT
80 INTEGER J, U11, INVD
81
82 COMPLEX*16 AK, AKKP1, AKP1, D, T
83 COMPLEX*16 U01_I_J, U01_IP1_J
84 COMPLEX*16 U11_I_J, U11_IP1_J
85 * ..
86 * .. External Functions ..
87 LOGICAL LSAME
88 EXTERNAL LSAME
89 * ..
90 * .. External Subroutines ..
91 EXTERNAL ZSYCONV, XERBLA, ZTRTRI
92 EXTERNAL ZGEMM, ZTRMM, ZHESWAPR
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC MAX
96 * ..
97 * .. Executable Statements ..
98 *
99 * Test the input parameters.
100 *
101 INFO = 0
102 UPPER = LSAME( UPLO, 'U' )
103 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
104 INFO = -1
105 ELSE IF( N.LT.0 ) THEN
106 INFO = -2
107 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
108 INFO = -4
109 END IF
110 *
111 * Quick return if possible
112 *
113 *
114 IF( INFO.NE.0 ) THEN
115 CALL XERBLA( 'ZHETRI2X', -INFO )
116 RETURN
117 END IF
118 IF( N.EQ.0 )
119 $ RETURN
120 *
121 * Convert A
122 * Workspace got Non-diag elements of D
123 *
124 CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
125 *
126 * Check that the diagonal matrix D is nonsingular.
127 *
128 IF( UPPER ) THEN
129 *
130 * Upper triangular storage: examine D from bottom to top
131 *
132 DO INFO = N, 1, -1
133 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
134 $ RETURN
135 END DO
136 ELSE
137 *
138 * Lower triangular storage: examine D from top to bottom.
139 *
140 DO INFO = 1, N
141 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
142 $ RETURN
143 END DO
144 END IF
145 INFO = 0
146 *
147 * Splitting Workspace
148 * U01 is a block (N,NB+1)
149 * The first element of U01 is in WORK(1,1)
150 * U11 is a block (NB+1,NB+1)
151 * The first element of U11 is in WORK(N+1,1)
152 U11 = N
153 * INVD is a block (N,2)
154 * The first element of INVD is in WORK(1,INVD)
155 INVD = NB+2
156
157 IF( UPPER ) THEN
158 *
159 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
160 *
161 CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
162 *
163 * inv(D) and inv(D)*inv(U)
164 *
165 K=1
166 DO WHILE ( K .LE. N )
167 IF( IPIV( K ).GT.0 ) THEN
168 * 1 x 1 diagonal NNB
169 WORK(K,INVD) = ONE / REAL ( A( K, K ) )
170 WORK(K,INVD+1) = 0
171 K=K+1
172 ELSE
173 * 2 x 2 diagonal NNB
174 T = ABS ( WORK(K+1,1) )
175 AK = REAL ( A( K, K ) ) / T
176 AKP1 = REAL ( A( K+1, K+1 ) ) / T
177 AKKP1 = WORK(K+1,1) / T
178 D = T*( AK*AKP1-ONE )
179 WORK(K,INVD) = AKP1 / D
180 WORK(K+1,INVD+1) = AK / D
181 WORK(K,INVD+1) = -AKKP1 / D
182 WORK(K+1,INVD) = DCONJG (WORK(K,INVD+1) )
183 K=K+2
184 END IF
185 END DO
186 *
187 * inv(U**H) = (inv(U))**H
188 *
189 * inv(U**H)*inv(D)*inv(U)
190 *
191 CUT=N
192 DO WHILE (CUT .GT. 0)
193 NNB=NB
194 IF (CUT .LE. NNB) THEN
195 NNB=CUT
196 ELSE
197 COUNT = 0
198 * count negative elements,
199 DO I=CUT+1-NNB,CUT
200 IF (IPIV(I) .LT. 0) COUNT=COUNT+1
201 END DO
202 * need a even number for a clear cut
203 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
204 END IF
205
206 CUT=CUT-NNB
207 *
208 * U01 Block
209 *
210 DO I=1,CUT
211 DO J=1,NNB
212 WORK(I,J)=A(I,CUT+J)
213 END DO
214 END DO
215 *
216 * U11 Block
217 *
218 DO I=1,NNB
219 WORK(U11+I,I)=CONE
220 DO J=1,I-1
221 WORK(U11+I,J)=ZERO
222 END DO
223 DO J=I+1,NNB
224 WORK(U11+I,J)=A(CUT+I,CUT+J)
225 END DO
226 END DO
227 *
228 * invD*U01
229 *
230 I=1
231 DO WHILE (I .LE. CUT)
232 IF (IPIV(I) > 0) THEN
233 DO J=1,NNB
234 WORK(I,J)=WORK(I,INVD)*WORK(I,J)
235 END DO
236 I=I+1
237 ELSE
238 DO J=1,NNB
239 U01_I_J = WORK(I,J)
240 U01_IP1_J = WORK(I+1,J)
241 WORK(I,J)=WORK(I,INVD)*U01_I_J+
242 $ WORK(I,INVD+1)*U01_IP1_J
243 WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
244 $ WORK(I+1,INVD+1)*U01_IP1_J
245 END DO
246 I=I+2
247 END IF
248 END DO
249 *
250 * invD1*U11
251 *
252 I=1
253 DO WHILE (I .LE. NNB)
254 IF (IPIV(CUT+I) > 0) THEN
255 DO J=I,NNB
256 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
257 END DO
258 I=I+1
259 ELSE
260 DO J=I,NNB
261 U11_I_J = WORK(U11+I,J)
262 U11_IP1_J = WORK(U11+I+1,J)
263 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
264 $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
265 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
266 $ WORK(CUT+I+1,INVD+1)*U11_IP1_J
267 END DO
268 I=I+2
269 END IF
270 END DO
271 *
272 * U11**H*invD1*U11->U11
273 *
274 CALL ZTRMM('L','U','C','U',NNB, NNB,
275 $ CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
276 *
277 DO I=1,NNB
278 DO J=I,NNB
279 A(CUT+I,CUT+J)=WORK(U11+I,J)
280 END DO
281 END DO
282 *
283 * U01**H*invD*U01->A(CUT+I,CUT+J)
284 *
285 CALL ZGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA,
286 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
287 *
288 * U11 = U11**H*invD1*U11 + U01**H*invD*U01
289 *
290 DO I=1,NNB
291 DO J=I,NNB
292 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
293 END DO
294 END DO
295 *
296 * U01 = U00**H*invD0*U01
297 *
298 CALL ZTRMM('L',UPLO,'C','U',CUT, NNB,
299 $ CONE,A,LDA,WORK,N+NB+1)
300
301 *
302 * Update U01
303 *
304 DO I=1,CUT
305 DO J=1,NNB
306 A(I,CUT+J)=WORK(I,J)
307 END DO
308 END DO
309 *
310 * Next Block
311 *
312 END DO
313 *
314 * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
315 *
316 I=1
317 DO WHILE ( I .LE. N )
318 IF( IPIV(I) .GT. 0 ) THEN
319 IP=IPIV(I)
320 IF (I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP )
321 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I )
322 ELSE
323 IP=-IPIV(I)
324 I=I+1
325 IF ( (I-1) .LT. IP)
326 $ CALL ZHESWAPR( UPLO, N, A, LDA, I-1 ,IP )
327 IF ( (I-1) .GT. IP)
328 $ CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I-1 )
329 ENDIF
330 I=I+1
331 END DO
332 ELSE
333 *
334 * LOWER...
335 *
336 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
337 *
338 CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
339 *
340 * inv(D) and inv(D)*inv(U)
341 *
342 K=N
343 DO WHILE ( K .GE. 1 )
344 IF( IPIV( K ).GT.0 ) THEN
345 * 1 x 1 diagonal NNB
346 WORK(K,INVD) = ONE / REAL ( A( K, K ) )
347 WORK(K,INVD+1) = 0
348 K=K-1
349 ELSE
350 * 2 x 2 diagonal NNB
351 T = ABS ( WORK(K-1,1) )
352 AK = REAL ( A( K-1, K-1 ) ) / T
353 AKP1 = REAL ( A( K, K ) ) / T
354 AKKP1 = WORK(K-1,1) / T
355 D = T*( AK*AKP1-ONE )
356 WORK(K-1,INVD) = AKP1 / D
357 WORK(K,INVD) = AK / D
358 WORK(K,INVD+1) = -AKKP1 / D
359 WORK(K-1,INVD+1) = DCONJG (WORK(K,INVD+1) )
360 K=K-2
361 END IF
362 END DO
363 *
364 * inv(U**H) = (inv(U))**H
365 *
366 * inv(U**H)*inv(D)*inv(U)
367 *
368 CUT=0
369 DO WHILE (CUT .LT. N)
370 NNB=NB
371 IF (CUT + NNB .GE. N) THEN
372 NNB=N-CUT
373 ELSE
374 COUNT = 0
375 * count negative elements,
376 DO I=CUT+1,CUT+NNB
377 IF (IPIV(I) .LT. 0) COUNT=COUNT+1
378 END DO
379 * need a even number for a clear cut
380 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
381 END IF
382 * L21 Block
383 DO I=1,N-CUT-NNB
384 DO J=1,NNB
385 WORK(I,J)=A(CUT+NNB+I,CUT+J)
386 END DO
387 END DO
388 * L11 Block
389 DO I=1,NNB
390 WORK(U11+I,I)=CONE
391 DO J=I+1,NNB
392 WORK(U11+I,J)=ZERO
393 END DO
394 DO J=1,I-1
395 WORK(U11+I,J)=A(CUT+I,CUT+J)
396 END DO
397 END DO
398 *
399 * invD*L21
400 *
401 I=N-CUT-NNB
402 DO WHILE (I .GE. 1)
403 IF (IPIV(CUT+NNB+I) > 0) THEN
404 DO J=1,NNB
405 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
406 END DO
407 I=I-1
408 ELSE
409 DO J=1,NNB
410 U01_I_J = WORK(I,J)
411 U01_IP1_J = WORK(I-1,J)
412 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
413 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
414 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
415 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
416 END DO
417 I=I-2
418 END IF
419 END DO
420 *
421 * invD1*L11
422 *
423 I=NNB
424 DO WHILE (I .GE. 1)
425 IF (IPIV(CUT+I) > 0) THEN
426 DO J=1,NNB
427 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
428 END DO
429 I=I-1
430 ELSE
431 DO J=1,NNB
432 U11_I_J = WORK(U11+I,J)
433 U11_IP1_J = WORK(U11+I-1,J)
434 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
435 $ WORK(CUT+I,INVD+1)*U11_IP1_J
436 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
437 $ WORK(CUT+I-1,INVD)*U11_IP1_J
438 END DO
439 I=I-2
440 END IF
441 END DO
442 *
443 * L11**H*invD1*L11->L11
444 *
445 CALL ZTRMM('L',UPLO,'C','U',NNB, NNB,
446 $ CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
447 *
448 DO I=1,NNB
449 DO J=1,I
450 A(CUT+I,CUT+J)=WORK(U11+I,J)
451 END DO
452 END DO
453 *
454 IF ( (CUT+NNB) .LT. N ) THEN
455 *
456 * L21**H*invD2*L21->A(CUT+I,CUT+J)
457 *
458 CALL ZGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1)
459 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
460
461 *
462 * L11 = L11**H*invD1*L11 + U01**H*invD*U01
463 *
464 DO I=1,NNB
465 DO J=1,I
466 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
467 END DO
468 END DO
469 *
470 * L01 = L22**H*invD2*L21
471 *
472 CALL ZTRMM('L',UPLO,'C','U', N-NNB-CUT, NNB,
473 $ CONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
474
475 * Update L21
476 DO I=1,N-CUT-NNB
477 DO J=1,NNB
478 A(CUT+NNB+I,CUT+J)=WORK(I,J)
479 END DO
480 END DO
481 ELSE
482 *
483 * L11 = L11**H*invD1*L11
484 *
485 DO I=1,NNB
486 DO J=1,I
487 A(CUT+I,CUT+J)=WORK(U11+I,J)
488 END DO
489 END DO
490 END IF
491 *
492 * Next Block
493 *
494 CUT=CUT+NNB
495 END DO
496 *
497 * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
498 *
499 I=N
500 DO WHILE ( I .GE. 1 )
501 IF( IPIV(I) .GT. 0 ) THEN
502 IP=IPIV(I)
503 IF (I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP )
504 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I )
505 ELSE
506 IP=-IPIV(I)
507 IF ( I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP )
508 IF ( I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I )
509 I=I-1
510 ENDIF
511 I=I-1
512 END DO
513 END IF
514 *
515 RETURN
516 *
517 * End of ZHETRI2X
518 *
519 END
520