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