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