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