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