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