1: SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
2: *
3: * -- LAPACK routine (version 3.3.0) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2010
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')*inv(D)*inv(U)*P'.
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') = (inv(U))'
186: *
187: * inv(U')*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: * U11T*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: * U01'invD*U01->A(CUT+I,CUT+J)
276: *
277: CALL ZGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
278: $ WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
279: *
280: * U11 = U11T*invD1*U11 + U01'invD*U01
281: *
282: DO I=1,NNB
283: DO J=I,NNB
284: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
285: END DO
286: END DO
287: *
288: * U01 = U00T*invD0*U01
289: *
290: CALL ZTRMM('L',UPLO,'T','U',CUT, NNB,
291: $ ONE,A,LDA,WORK,N+NB+1)
292:
293: *
294: * Update U01
295: *
296: DO I=1,CUT
297: DO J=1,NNB
298: A(I,CUT+J)=WORK(I,J)
299: END DO
300: END DO
301: *
302: * Next Block
303: *
304: END DO
305: *
306: * Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
307: *
308: I=1
309: DO WHILE ( I .LE. N )
310: IF( IPIV(I) .GT. 0 ) THEN
311: IP=IPIV(I)
312: IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP )
313: IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )
314: ELSE
315: IP=-IPIV(I)
316: I=I+1
317: IF ( (I-1) .LT. IP)
318: $ CALL ZSYSWAPR( UPLO, N, A, I-1 ,IP )
319: IF ( (I-1) .GT. IP)
320: $ CALL ZSYSWAPR( UPLO, N, A, IP ,I-1 )
321: ENDIF
322: I=I+1
323: END DO
324: ELSE
325: *
326: * LOWER...
327: *
328: * invA = P * inv(U')*inv(D)*inv(U)*P'.
329: *
330: CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
331: *
332: * inv(D) and inv(D)*inv(U)
333: *
334: K=N
335: DO WHILE ( K .GE. 1 )
336: IF( IPIV( K ).GT.0 ) THEN
337: * 1 x 1 diagonal NNB
338: WORK(K,INVD) = 1/ A( K, K )
339: WORK(K,INVD+1) = 0
340: K=K-1
341: ELSE
342: * 2 x 2 diagonal NNB
343: T = WORK(K-1,1)
344: AK = A( K-1, K-1 ) / T
345: AKP1 = A( K, K ) / T
346: AKKP1 = WORK(K-1,1) / T
347: D = T*( AK*AKP1-ONE )
348: WORK(K-1,INVD) = AKP1 / D
349: WORK(K,INVD) = AK / D
350: WORK(K,INVD+1) = -AKKP1 / D
351: WORK(K-1,INVD+1) = -AKKP1 / D
352: K=K-2
353: END IF
354: END DO
355: *
356: * inv(U') = (inv(U))'
357: *
358: * inv(U')*inv(D)*inv(U)
359: *
360: CUT=0
361: DO WHILE (CUT .LT. N)
362: NNB=NB
363: IF (CUT + NNB .GE. N) THEN
364: NNB=N-CUT
365: ELSE
366: COUNT = 0
367: * count negative elements,
368: DO I=CUT+1,CUT+NNB
369: IF (IPIV(I) .LT. 0) COUNT=COUNT+1
370: END DO
371: * need a even number for a clear cut
372: IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
373: END IF
374: * L21 Block
375: DO I=1,N-CUT-NNB
376: DO J=1,NNB
377: WORK(I,J)=A(CUT+NNB+I,CUT+J)
378: END DO
379: END DO
380: * L11 Block
381: DO I=1,NNB
382: WORK(U11+I,I)=ONE
383: DO J=I+1,NNB
384: WORK(U11+I,J)=ZERO
385: END DO
386: DO J=1,I-1
387: WORK(U11+I,J)=A(CUT+I,CUT+J)
388: END DO
389: END DO
390: *
391: * invD*L21
392: *
393: I=N-CUT-NNB
394: DO WHILE (I .GE. 1)
395: IF (IPIV(CUT+NNB+I) > 0) THEN
396: DO J=1,NNB
397: WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
398: END DO
399: I=I-1
400: ELSE
401: DO J=1,NNB
402: U01_I_J = WORK(I,J)
403: U01_IP1_J = WORK(I-1,J)
404: WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
405: $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
406: WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
407: $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
408: END DO
409: I=I-2
410: END IF
411: END DO
412: *
413: * invD1*L11
414: *
415: I=NNB
416: DO WHILE (I .GE. 1)
417: IF (IPIV(CUT+I) > 0) THEN
418: DO J=1,NNB
419: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
420: END DO
421: I=I-1
422: ELSE
423: DO J=1,NNB
424: U11_I_J = WORK(U11+I,J)
425: U11_IP1_J = WORK(U11+I-1,J)
426: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
427: $ WORK(CUT+I,INVD+1)*U11_IP1_J
428: WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
429: $ WORK(CUT+I-1,INVD)*U11_IP1_J
430: END DO
431: I=I-2
432: END IF
433: END DO
434: *
435: * L11T*invD1*L11->L11
436: *
437: CALL ZTRMM('L',UPLO,'T','U',NNB, NNB,
438: $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
439:
440: IF ( (CUT+NNB) .LT. N ) THEN
441: *
442: * L21T*invD2*L21->A(CUT+I,CUT+J)
443: *
444: CALL ZGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
445: $ ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
446:
447: *
448: * L11 = L11T*invD1*L11 + U01'invD*U01
449: *
450: DO I=1,NNB
451: DO J=1,I
452: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
453: END DO
454: END DO
455: *
456: * U01 = L22T*invD2*L21
457: *
458: CALL ZTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
459: $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
460:
461: * Update L21
462: DO I=1,N-CUT-NNB
463: DO J=1,NNB
464: A(CUT+NNB+I,CUT+J)=WORK(I,J)
465: END DO
466: END DO
467: ELSE
468: *
469: * L11 = L11T*invD1*L11
470: *
471: DO I=1,NNB
472: DO J=1,I
473: A(CUT+I,CUT+J)=WORK(U11+I,J)
474: END DO
475: END DO
476: END IF
477: *
478: * Next Block
479: *
480: CUT=CUT+NNB
481: END DO
482: *
483: * Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
484: *
485: I=N
486: DO WHILE ( I .GE. 1 )
487: IF( IPIV(I) .GT. 0 ) THEN
488: IP=IPIV(I)
489: IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP )
490: IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )
491: ELSE
492: IP=-IPIV(I)
493: IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP )
494: IF ( I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )
495: I=I-1
496: ENDIF
497: I=I-1
498: END DO
499: END IF
500: *
501: RETURN
502: *
503: * End of ZSYTRI2X
504: *
505: END
506:
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