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:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetri2x.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack