1: *> \brief \b ZHETRI_ROOK computes the inverse of HE matrix using the factorization obtained with the bounded Bunch-Kaufman ("rook") diagonal pivoting method.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZHETRI_ROOK + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetri_rook.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetri_rook.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetri_rook.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHETRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, LDA, N
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * COMPLEX*16 A( LDA, * ), WORK( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZHETRI_ROOK computes the inverse of a complex Hermitian indefinite matrix
39: *> A using the factorization A = U*D*U**H or A = L*D*L**H computed by
40: *> ZHETRF_ROOK.
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 block diagonal matrix D and the multipliers
65: *> used to obtain the factor U or L as computed by ZHETRF_ROOK.
66: *>
67: *> On exit, if INFO = 0, the (Hermitian) 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 block structure of D
85: *> as determined by ZHETRF_ROOK.
86: *> \endverbatim
87: *>
88: *> \param[out] WORK
89: *> \verbatim
90: *> WORK is COMPLEX*16 array, dimension (N)
91: *> \endverbatim
92: *>
93: *> \param[out] INFO
94: *> \verbatim
95: *> INFO is INTEGER
96: *> = 0: successful exit
97: *> < 0: if INFO = -i, the i-th argument had an illegal value
98: *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
99: *> inverse could not be computed.
100: *> \endverbatim
101: *
102: * Authors:
103: * ========
104: *
105: *> \author Univ. of Tennessee
106: *> \author Univ. of California Berkeley
107: *> \author Univ. of Colorado Denver
108: *> \author NAG Ltd.
109: *
110: *> \ingroup complex16HEcomputational
111: *
112: *> \par Contributors:
113: * ==================
114: *>
115: *> \verbatim
116: *>
117: *> November 2013, Igor Kozachenko,
118: *> Computer Science Division,
119: *> University of California, Berkeley
120: *>
121: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
122: *> School of Mathematics,
123: *> University of Manchester
124: *> \endverbatim
125: *
126: * =====================================================================
127: SUBROUTINE ZHETRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
128: *
129: * -- LAPACK computational routine --
130: * -- LAPACK is a software package provided by Univ. of Tennessee, --
131: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
132: *
133: * .. Scalar Arguments ..
134: CHARACTER UPLO
135: INTEGER INFO, LDA, N
136: * ..
137: * .. Array Arguments ..
138: INTEGER IPIV( * )
139: COMPLEX*16 A( LDA, * ), WORK( * )
140: * ..
141: *
142: * =====================================================================
143: *
144: * .. Parameters ..
145: DOUBLE PRECISION ONE
146: COMPLEX*16 CONE, CZERO
147: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
148: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
149: * ..
150: * .. Local Scalars ..
151: LOGICAL UPPER
152: INTEGER J, K, KP, KSTEP
153: DOUBLE PRECISION AK, AKP1, D, T
154: COMPLEX*16 AKKP1, TEMP
155: * ..
156: * .. External Functions ..
157: LOGICAL LSAME
158: COMPLEX*16 ZDOTC
159: EXTERNAL LSAME, ZDOTC
160: * ..
161: * .. External Subroutines ..
162: EXTERNAL ZCOPY, ZHEMV, ZSWAP, XERBLA
163: * ..
164: * .. Intrinsic Functions ..
165: INTRINSIC ABS, DCONJG, MAX, DBLE
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: IF( INFO.NE.0 ) THEN
181: CALL XERBLA( 'ZHETRI_ROOK', -INFO )
182: RETURN
183: END IF
184: *
185: * Quick return if possible
186: *
187: IF( N.EQ.0 )
188: $ RETURN
189: *
190: * Check that the diagonal matrix D is nonsingular.
191: *
192: IF( UPPER ) THEN
193: *
194: * Upper triangular storage: examine D from bottom to top
195: *
196: DO 10 INFO = N, 1, -1
197: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
198: $ RETURN
199: 10 CONTINUE
200: ELSE
201: *
202: * Lower triangular storage: examine D from top to bottom.
203: *
204: DO 20 INFO = 1, N
205: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
206: $ RETURN
207: 20 CONTINUE
208: END IF
209: INFO = 0
210: *
211: IF( UPPER ) THEN
212: *
213: * Compute inv(A) from the factorization A = U*D*U**H.
214: *
215: * K is the main loop index, increasing from 1 to N in steps of
216: * 1 or 2, depending on the size of the diagonal blocks.
217: *
218: K = 1
219: 30 CONTINUE
220: *
221: * If K > N, exit from loop.
222: *
223: IF( K.GT.N )
224: $ GO TO 70
225: *
226: IF( IPIV( K ).GT.0 ) THEN
227: *
228: * 1 x 1 diagonal block
229: *
230: * Invert the diagonal block.
231: *
232: A( K, K ) = ONE / DBLE( A( K, K ) )
233: *
234: * Compute column K of the inverse.
235: *
236: IF( K.GT.1 ) THEN
237: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
238: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
239: $ A( 1, K ), 1 )
240: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
241: $ K ), 1 ) )
242: END IF
243: KSTEP = 1
244: ELSE
245: *
246: * 2 x 2 diagonal block
247: *
248: * Invert the diagonal block.
249: *
250: T = ABS( A( K, K+1 ) )
251: AK = DBLE( A( K, K ) ) / T
252: AKP1 = DBLE( A( K+1, K+1 ) ) / T
253: AKKP1 = A( K, K+1 ) / T
254: D = T*( AK*AKP1-ONE )
255: A( K, K ) = AKP1 / D
256: A( K+1, K+1 ) = AK / D
257: A( K, K+1 ) = -AKKP1 / D
258: *
259: * Compute columns K and K+1 of the inverse.
260: *
261: IF( K.GT.1 ) THEN
262: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
263: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
264: $ A( 1, K ), 1 )
265: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
266: $ K ), 1 ) )
267: A( K, K+1 ) = A( K, K+1 ) -
268: $ ZDOTC( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
269: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
270: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
271: $ A( 1, K+1 ), 1 )
272: A( K+1, K+1 ) = A( K+1, K+1 ) -
273: $ DBLE( ZDOTC( K-1, WORK, 1, A( 1, K+1 ),
274: $ 1 ) )
275: END IF
276: KSTEP = 2
277: END IF
278: *
279: IF( KSTEP.EQ.1 ) THEN
280: *
281: * Interchange rows and columns K and IPIV(K) in the leading
282: * submatrix A(1:k,1:k)
283: *
284: KP = IPIV( K )
285: IF( KP.NE.K ) THEN
286: *
287: IF( KP.GT.1 )
288: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
289: *
290: DO 40 J = KP + 1, K - 1
291: TEMP = DCONJG( A( J, K ) )
292: A( J, K ) = DCONJG( A( KP, J ) )
293: A( KP, J ) = TEMP
294: 40 CONTINUE
295: *
296: A( KP, K ) = DCONJG( A( KP, K ) )
297: *
298: TEMP = A( K, K )
299: A( K, K ) = A( KP, KP )
300: A( KP, KP ) = TEMP
301: END IF
302: ELSE
303: *
304: * Interchange rows and columns K and K+1 with -IPIV(K) and
305: * -IPIV(K+1) in the leading submatrix A(k+1:n,k+1:n)
306: *
307: * (1) Interchange rows and columns K and -IPIV(K)
308: *
309: KP = -IPIV( K )
310: IF( KP.NE.K ) THEN
311: *
312: IF( KP.GT.1 )
313: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
314: *
315: DO 50 J = KP + 1, K - 1
316: TEMP = DCONJG( A( J, K ) )
317: A( J, K ) = DCONJG( A( KP, J ) )
318: A( KP, J ) = TEMP
319: 50 CONTINUE
320: *
321: A( KP, K ) = DCONJG( A( KP, K ) )
322: *
323: TEMP = A( K, K )
324: A( K, K ) = A( KP, KP )
325: A( KP, KP ) = TEMP
326: *
327: TEMP = A( K, K+1 )
328: A( K, K+1 ) = A( KP, K+1 )
329: A( KP, K+1 ) = TEMP
330: END IF
331: *
332: * (2) Interchange rows and columns K+1 and -IPIV(K+1)
333: *
334: K = K + 1
335: KP = -IPIV( K )
336: IF( KP.NE.K ) THEN
337: *
338: IF( KP.GT.1 )
339: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
340: *
341: DO 60 J = KP + 1, K - 1
342: TEMP = DCONJG( A( J, K ) )
343: A( J, K ) = DCONJG( A( KP, J ) )
344: A( KP, J ) = TEMP
345: 60 CONTINUE
346: *
347: A( KP, K ) = DCONJG( A( KP, K ) )
348: *
349: TEMP = A( K, K )
350: A( K, K ) = A( KP, KP )
351: A( KP, KP ) = TEMP
352: END IF
353: END IF
354: *
355: K = K + 1
356: GO TO 30
357: 70 CONTINUE
358: *
359: ELSE
360: *
361: * Compute inv(A) from the factorization A = L*D*L**H.
362: *
363: * K is the main loop index, decreasing from N to 1 in steps of
364: * 1 or 2, depending on the size of the diagonal blocks.
365: *
366: K = N
367: 80 CONTINUE
368: *
369: * If K < 1, exit from loop.
370: *
371: IF( K.LT.1 )
372: $ GO TO 120
373: *
374: IF( IPIV( K ).GT.0 ) THEN
375: *
376: * 1 x 1 diagonal block
377: *
378: * Invert the diagonal block.
379: *
380: A( K, K ) = ONE / DBLE( A( K, K ) )
381: *
382: * Compute column K of the inverse.
383: *
384: IF( K.LT.N ) THEN
385: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
386: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
387: $ 1, CZERO, A( K+1, K ), 1 )
388: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
389: $ A( K+1, K ), 1 ) )
390: END IF
391: KSTEP = 1
392: ELSE
393: *
394: * 2 x 2 diagonal block
395: *
396: * Invert the diagonal block.
397: *
398: T = ABS( A( K, K-1 ) )
399: AK = DBLE( A( K-1, K-1 ) ) / T
400: AKP1 = DBLE( A( K, K ) ) / T
401: AKKP1 = A( K, K-1 ) / T
402: D = T*( AK*AKP1-ONE )
403: A( K-1, K-1 ) = AKP1 / D
404: A( K, K ) = AK / D
405: A( K, K-1 ) = -AKKP1 / D
406: *
407: * Compute columns K-1 and K of the inverse.
408: *
409: IF( K.LT.N ) THEN
410: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
411: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
412: $ 1, CZERO, A( K+1, K ), 1 )
413: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
414: $ A( K+1, K ), 1 ) )
415: A( K, K-1 ) = A( K, K-1 ) -
416: $ ZDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
417: $ 1 )
418: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
419: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
420: $ 1, CZERO, A( K+1, K-1 ), 1 )
421: A( K-1, K-1 ) = A( K-1, K-1 ) -
422: $ DBLE( ZDOTC( N-K, WORK, 1, A( K+1, K-1 ),
423: $ 1 ) )
424: END IF
425: KSTEP = 2
426: END IF
427: *
428: IF( KSTEP.EQ.1 ) THEN
429: *
430: * Interchange rows and columns K and IPIV(K) in the trailing
431: * submatrix A(k:n,k:n)
432: *
433: KP = IPIV( K )
434: IF( KP.NE.K ) THEN
435: *
436: IF( KP.LT.N )
437: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
438: *
439: DO 90 J = K + 1, KP - 1
440: TEMP = DCONJG( A( J, K ) )
441: A( J, K ) = DCONJG( A( KP, J ) )
442: A( KP, J ) = TEMP
443: 90 CONTINUE
444: *
445: A( KP, K ) = DCONJG( A( KP, K ) )
446: *
447: TEMP = A( K, K )
448: A( K, K ) = A( KP, KP )
449: A( KP, KP ) = TEMP
450: END IF
451: ELSE
452: *
453: * Interchange rows and columns K and K-1 with -IPIV(K) and
454: * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
455: *
456: * (1) Interchange rows and columns K and -IPIV(K)
457: *
458: KP = -IPIV( K )
459: IF( KP.NE.K ) THEN
460: *
461: IF( KP.LT.N )
462: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
463: *
464: DO 100 J = K + 1, KP - 1
465: TEMP = DCONJG( A( J, K ) )
466: A( J, K ) = DCONJG( A( KP, J ) )
467: A( KP, J ) = TEMP
468: 100 CONTINUE
469: *
470: A( KP, K ) = DCONJG( A( KP, K ) )
471: *
472: TEMP = A( K, K )
473: A( K, K ) = A( KP, KP )
474: A( KP, KP ) = TEMP
475: *
476: TEMP = A( K, K-1 )
477: A( K, K-1 ) = A( KP, K-1 )
478: A( KP, K-1 ) = TEMP
479: END IF
480: *
481: * (2) Interchange rows and columns K-1 and -IPIV(K-1)
482: *
483: K = K - 1
484: KP = -IPIV( K )
485: IF( KP.NE.K ) THEN
486: *
487: IF( KP.LT.N )
488: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
489: *
490: DO 110 J = K + 1, KP - 1
491: TEMP = DCONJG( A( J, K ) )
492: A( J, K ) = DCONJG( A( KP, J ) )
493: A( KP, J ) = TEMP
494: 110 CONTINUE
495: *
496: A( KP, K ) = DCONJG( A( KP, K ) )
497: *
498: TEMP = A( K, K )
499: A( K, K ) = A( KP, KP )
500: A( KP, KP ) = TEMP
501: END IF
502: END IF
503: *
504: K = K - 1
505: GO TO 80
506: 120 CONTINUE
507: END IF
508: *
509: RETURN
510: *
511: * End of ZHETRI_ROOK
512: *
513: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetri_rook.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack