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