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+NB+1,NB+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: *> \ingroup complex16SYcomputational
117: *
118: * =====================================================================
119: SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
120: *
121: * -- LAPACK computational routine --
122: * -- LAPACK is a software package provided by Univ. of Tennessee, --
123: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124: *
125: * .. Scalar Arguments ..
126: CHARACTER UPLO
127: INTEGER INFO, LDA, N, NB
128: * ..
129: * .. Array Arguments ..
130: INTEGER IPIV( * )
131: COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* )
132: * ..
133: *
134: * =====================================================================
135: *
136: * .. Parameters ..
137: COMPLEX*16 ONE, ZERO
138: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
139: $ ZERO = ( 0.0D+0, 0.0D+0 ) )
140: * ..
141: * .. Local Scalars ..
142: LOGICAL UPPER
143: INTEGER I, IINFO, IP, K, CUT, NNB
144: INTEGER COUNT
145: INTEGER J, U11, INVD
146:
147: COMPLEX*16 AK, AKKP1, AKP1, D, T
148: COMPLEX*16 U01_I_J, U01_IP1_J
149: COMPLEX*16 U11_I_J, U11_IP1_J
150: * ..
151: * .. External Functions ..
152: LOGICAL LSAME
153: EXTERNAL LSAME
154: * ..
155: * .. External Subroutines ..
156: EXTERNAL ZSYCONV, XERBLA, ZTRTRI
157: EXTERNAL ZGEMM, ZTRMM, ZSYSWAPR
158: * ..
159: * .. Intrinsic Functions ..
160: INTRINSIC MAX
161: * ..
162: * .. Executable Statements ..
163: *
164: * Test the input parameters.
165: *
166: INFO = 0
167: UPPER = LSAME( UPLO, 'U' )
168: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
169: INFO = -1
170: ELSE IF( N.LT.0 ) THEN
171: INFO = -2
172: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
173: INFO = -4
174: END IF
175: *
176: * Quick return if possible
177: *
178: *
179: IF( INFO.NE.0 ) THEN
180: CALL XERBLA( 'ZSYTRI2X', -INFO )
181: RETURN
182: END IF
183: IF( N.EQ.0 )
184: $ RETURN
185: *
186: * Convert A
187: * Workspace got Non-diag elements of D
188: *
189: CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
190: *
191: * Check that the diagonal matrix D is nonsingular.
192: *
193: IF( UPPER ) THEN
194: *
195: * Upper triangular storage: examine D from bottom to top
196: *
197: DO INFO = N, 1, -1
198: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
199: $ RETURN
200: END DO
201: ELSE
202: *
203: * Lower triangular storage: examine D from top to bottom.
204: *
205: DO INFO = 1, N
206: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
207: $ RETURN
208: END DO
209: END IF
210: INFO = 0
211: *
212: * Splitting Workspace
213: * U01 is a block (N,NB+1)
214: * The first element of U01 is in WORK(1,1)
215: * U11 is a block (NB+1,NB+1)
216: * The first element of U11 is in WORK(N+1,1)
217: U11 = N
218: * INVD is a block (N,2)
219: * The first element of INVD is in WORK(1,INVD)
220: INVD = NB+2
221:
222: IF( UPPER ) THEN
223: *
224: * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
225: *
226: CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
227: *
228: * inv(D) and inv(D)*inv(U)
229: *
230: K=1
231: DO WHILE ( K .LE. N )
232: IF( IPIV( K ).GT.0 ) THEN
233: * 1 x 1 diagonal NNB
234: WORK(K,INVD) = ONE / A( K, K )
235: WORK(K,INVD+1) = 0
236: K=K+1
237: ELSE
238: * 2 x 2 diagonal NNB
239: T = WORK(K+1,1)
240: AK = A( K, K ) / T
241: AKP1 = A( K+1, K+1 ) / T
242: AKKP1 = WORK(K+1,1) / T
243: D = T*( AK*AKP1-ONE )
244: WORK(K,INVD) = AKP1 / D
245: WORK(K+1,INVD+1) = AK / D
246: WORK(K,INVD+1) = -AKKP1 / D
247: WORK(K+1,INVD) = -AKKP1 / D
248: K=K+2
249: END IF
250: END DO
251: *
252: * inv(U**T) = (inv(U))**T
253: *
254: * inv(U**T)*inv(D)*inv(U)
255: *
256: CUT=N
257: DO WHILE (CUT .GT. 0)
258: NNB=NB
259: IF (CUT .LE. NNB) THEN
260: NNB=CUT
261: ELSE
262: COUNT = 0
263: * count negative elements,
264: DO I=CUT+1-NNB,CUT
265: IF (IPIV(I) .LT. 0) COUNT=COUNT+1
266: END DO
267: * need a even number for a clear cut
268: IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
269: END IF
270:
271: CUT=CUT-NNB
272: *
273: * U01 Block
274: *
275: DO I=1,CUT
276: DO J=1,NNB
277: WORK(I,J)=A(I,CUT+J)
278: END DO
279: END DO
280: *
281: * U11 Block
282: *
283: DO I=1,NNB
284: WORK(U11+I,I)=ONE
285: DO J=1,I-1
286: WORK(U11+I,J)=ZERO
287: END DO
288: DO J=I+1,NNB
289: WORK(U11+I,J)=A(CUT+I,CUT+J)
290: END DO
291: END DO
292: *
293: * invD*U01
294: *
295: I=1
296: DO WHILE (I .LE. CUT)
297: IF (IPIV(I) > 0) THEN
298: DO J=1,NNB
299: WORK(I,J)=WORK(I,INVD)*WORK(I,J)
300: END DO
301: I=I+1
302: ELSE
303: DO J=1,NNB
304: U01_I_J = WORK(I,J)
305: U01_IP1_J = WORK(I+1,J)
306: WORK(I,J)=WORK(I,INVD)*U01_I_J+
307: $ WORK(I,INVD+1)*U01_IP1_J
308: WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
309: $ WORK(I+1,INVD+1)*U01_IP1_J
310: END DO
311: I=I+2
312: END IF
313: END DO
314: *
315: * invD1*U11
316: *
317: I=1
318: DO WHILE (I .LE. NNB)
319: IF (IPIV(CUT+I) > 0) THEN
320: DO J=I,NNB
321: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
322: END DO
323: I=I+1
324: ELSE
325: DO J=I,NNB
326: U11_I_J = WORK(U11+I,J)
327: U11_IP1_J = WORK(U11+I+1,J)
328: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
329: $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
330: WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
331: $ WORK(CUT+I+1,INVD+1)*U11_IP1_J
332: END DO
333: I=I+2
334: END IF
335: END DO
336: *
337: * U11**T*invD1*U11->U11
338: *
339: CALL ZTRMM('L','U','T','U',NNB, NNB,
340: $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
341: *
342: DO I=1,NNB
343: DO J=I,NNB
344: A(CUT+I,CUT+J)=WORK(U11+I,J)
345: END DO
346: END DO
347: *
348: * U01**T*invD*U01->A(CUT+I,CUT+J)
349: *
350: CALL ZGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
351: $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
352: *
353: * U11 = U11**T*invD1*U11 + U01**T*invD*U01
354: *
355: DO I=1,NNB
356: DO J=I,NNB
357: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
358: END DO
359: END DO
360: *
361: * U01 = U00**T*invD0*U01
362: *
363: CALL ZTRMM('L',UPLO,'T','U',CUT, NNB,
364: $ ONE,A,LDA,WORK,N+NB+1)
365:
366: *
367: * Update U01
368: *
369: DO I=1,CUT
370: DO J=1,NNB
371: A(I,CUT+J)=WORK(I,J)
372: END DO
373: END DO
374: *
375: * Next Block
376: *
377: END DO
378: *
379: * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
380: *
381: I=1
382: DO WHILE ( I .LE. N )
383: IF( IPIV(I) .GT. 0 ) THEN
384: IP=IPIV(I)
385: IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP )
386: IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
387: ELSE
388: IP=-IPIV(I)
389: I=I+1
390: IF ( (I-1) .LT. IP)
391: $ CALL ZSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
392: IF ( (I-1) .GT. IP)
393: $ CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I-1 )
394: ENDIF
395: I=I+1
396: END DO
397: ELSE
398: *
399: * LOWER...
400: *
401: * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
402: *
403: CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
404: *
405: * inv(D) and inv(D)*inv(U)
406: *
407: K=N
408: DO WHILE ( K .GE. 1 )
409: IF( IPIV( K ).GT.0 ) THEN
410: * 1 x 1 diagonal NNB
411: WORK(K,INVD) = ONE / A( K, K )
412: WORK(K,INVD+1) = 0
413: K=K-1
414: ELSE
415: * 2 x 2 diagonal NNB
416: T = WORK(K-1,1)
417: AK = A( K-1, K-1 ) / T
418: AKP1 = A( K, K ) / T
419: AKKP1 = WORK(K-1,1) / T
420: D = T*( AK*AKP1-ONE )
421: WORK(K-1,INVD) = AKP1 / D
422: WORK(K,INVD) = AK / D
423: WORK(K,INVD+1) = -AKKP1 / D
424: WORK(K-1,INVD+1) = -AKKP1 / D
425: K=K-2
426: END IF
427: END DO
428: *
429: * inv(U**T) = (inv(U))**T
430: *
431: * inv(U**T)*inv(D)*inv(U)
432: *
433: CUT=0
434: DO WHILE (CUT .LT. N)
435: NNB=NB
436: IF (CUT + NNB .GE. N) THEN
437: NNB=N-CUT
438: ELSE
439: COUNT = 0
440: * count negative elements,
441: DO I=CUT+1,CUT+NNB
442: IF (IPIV(I) .LT. 0) COUNT=COUNT+1
443: END DO
444: * need a even number for a clear cut
445: IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
446: END IF
447: * L21 Block
448: DO I=1,N-CUT-NNB
449: DO J=1,NNB
450: WORK(I,J)=A(CUT+NNB+I,CUT+J)
451: END DO
452: END DO
453: * L11 Block
454: DO I=1,NNB
455: WORK(U11+I,I)=ONE
456: DO J=I+1,NNB
457: WORK(U11+I,J)=ZERO
458: END DO
459: DO J=1,I-1
460: WORK(U11+I,J)=A(CUT+I,CUT+J)
461: END DO
462: END DO
463: *
464: * invD*L21
465: *
466: I=N-CUT-NNB
467: DO WHILE (I .GE. 1)
468: IF (IPIV(CUT+NNB+I) > 0) THEN
469: DO J=1,NNB
470: WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
471: END DO
472: I=I-1
473: ELSE
474: DO J=1,NNB
475: U01_I_J = WORK(I,J)
476: U01_IP1_J = WORK(I-1,J)
477: WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
478: $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
479: WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
480: $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
481: END DO
482: I=I-2
483: END IF
484: END DO
485: *
486: * invD1*L11
487: *
488: I=NNB
489: DO WHILE (I .GE. 1)
490: IF (IPIV(CUT+I) > 0) THEN
491: DO J=1,NNB
492: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
493: END DO
494: I=I-1
495: ELSE
496: DO J=1,NNB
497: U11_I_J = WORK(U11+I,J)
498: U11_IP1_J = WORK(U11+I-1,J)
499: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
500: $ WORK(CUT+I,INVD+1)*U11_IP1_J
501: WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
502: $ WORK(CUT+I-1,INVD)*U11_IP1_J
503: END DO
504: I=I-2
505: END IF
506: END DO
507: *
508: * L11**T*invD1*L11->L11
509: *
510: CALL ZTRMM('L',UPLO,'T','U',NNB, NNB,
511: $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
512: *
513: DO I=1,NNB
514: DO J=1,I
515: A(CUT+I,CUT+J)=WORK(U11+I,J)
516: END DO
517: END DO
518: *
519:
520: IF ( (CUT+NNB) .LT. N ) THEN
521: *
522: * L21**T*invD2*L21->A(CUT+I,CUT+J)
523: *
524: CALL ZGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
525: $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
526:
527: *
528: * L11 = L11**T*invD1*L11 + U01**T*invD*U01
529: *
530: DO I=1,NNB
531: DO J=1,I
532: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
533: END DO
534: END DO
535: *
536: * U01 = L22**T*invD2*L21
537: *
538: CALL ZTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
539: $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
540:
541: * Update L21
542: DO I=1,N-CUT-NNB
543: DO J=1,NNB
544: A(CUT+NNB+I,CUT+J)=WORK(I,J)
545: END DO
546: END DO
547: ELSE
548: *
549: * L11 = L11**T*invD1*L11
550: *
551: DO I=1,NNB
552: DO J=1,I
553: A(CUT+I,CUT+J)=WORK(U11+I,J)
554: END DO
555: END DO
556: END IF
557: *
558: * Next Block
559: *
560: CUT=CUT+NNB
561: END DO
562: *
563: * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
564: *
565: I=N
566: DO WHILE ( I .GE. 1 )
567: IF( IPIV(I) .GT. 0 ) THEN
568: IP=IPIV(I)
569: IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP )
570: IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
571: ELSE
572: IP=-IPIV(I)
573: IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP )
574: IF ( I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
575: I=I-1
576: ENDIF
577: I=I-1
578: END DO
579: END IF
580: *
581: RETURN
582: *
583: * End of ZSYTRI2X
584: *
585: END
586:
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