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: *> \date November 2013
111: *
112: *> \ingroup complex16HEcomputational
113: *
114: *> \par Contributors:
115: * ==================
116: *>
117: *> \verbatim
118: *>
119: *> November 2013, Igor Kozachenko,
120: *> Computer Science Division,
121: *> University of California, Berkeley
122: *>
123: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
124: *> School of Mathematics,
125: *> University of Manchester
126: *> \endverbatim
127: *
128: * =====================================================================
129: SUBROUTINE ZHETRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
130: *
131: * -- LAPACK computational routine (version 3.5.0) --
132: * -- LAPACK is a software package provided by Univ. of Tennessee, --
133: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134: * November 2013
135: *
136: * .. Scalar Arguments ..
137: CHARACTER UPLO
138: INTEGER INFO, LDA, N
139: * ..
140: * .. Array Arguments ..
141: INTEGER IPIV( * )
142: COMPLEX*16 A( LDA, * ), WORK( * )
143: * ..
144: *
145: * =====================================================================
146: *
147: * .. Parameters ..
148: DOUBLE PRECISION ONE
149: COMPLEX*16 CONE, CZERO
150: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
151: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
152: * ..
153: * .. Local Scalars ..
154: LOGICAL UPPER
155: INTEGER J, K, KP, KSTEP
156: DOUBLE PRECISION AK, AKP1, D, T
157: COMPLEX*16 AKKP1, TEMP
158: * ..
159: * .. External Functions ..
160: LOGICAL LSAME
161: COMPLEX*16 ZDOTC
162: EXTERNAL LSAME, ZDOTC
163: * ..
164: * .. External Subroutines ..
165: EXTERNAL ZCOPY, ZHEMV, ZSWAP, XERBLA
166: * ..
167: * .. Intrinsic Functions ..
168: INTRINSIC ABS, DCONJG, MAX, DBLE
169: * ..
170: * .. Executable Statements ..
171: *
172: * Test the input parameters.
173: *
174: INFO = 0
175: UPPER = LSAME( UPLO, 'U' )
176: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
177: INFO = -1
178: ELSE IF( N.LT.0 ) THEN
179: INFO = -2
180: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
181: INFO = -4
182: END IF
183: IF( INFO.NE.0 ) THEN
184: CALL XERBLA( 'ZHETRI_ROOK', -INFO )
185: RETURN
186: END IF
187: *
188: * Quick return if possible
189: *
190: IF( N.EQ.0 )
191: $ RETURN
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 10 INFO = N, 1, -1
200: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
201: $ RETURN
202: 10 CONTINUE
203: ELSE
204: *
205: * Lower triangular storage: examine D from top to bottom.
206: *
207: DO 20 INFO = 1, N
208: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
209: $ RETURN
210: 20 CONTINUE
211: END IF
212: INFO = 0
213: *
214: IF( UPPER ) THEN
215: *
216: * Compute inv(A) from the factorization A = U*D*U**H.
217: *
218: * K is the main loop index, increasing from 1 to N in steps of
219: * 1 or 2, depending on the size of the diagonal blocks.
220: *
221: K = 1
222: 30 CONTINUE
223: *
224: * If K > N, exit from loop.
225: *
226: IF( K.GT.N )
227: $ GO TO 70
228: *
229: IF( IPIV( K ).GT.0 ) THEN
230: *
231: * 1 x 1 diagonal block
232: *
233: * Invert the diagonal block.
234: *
235: A( K, K ) = ONE / DBLE( A( K, K ) )
236: *
237: * Compute column K of the inverse.
238: *
239: IF( K.GT.1 ) THEN
240: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
241: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
242: $ A( 1, K ), 1 )
243: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
244: $ K ), 1 ) )
245: END IF
246: KSTEP = 1
247: ELSE
248: *
249: * 2 x 2 diagonal block
250: *
251: * Invert the diagonal block.
252: *
253: T = ABS( A( K, K+1 ) )
254: AK = DBLE( A( K, K ) ) / T
255: AKP1 = DBLE( A( K+1, K+1 ) ) / T
256: AKKP1 = A( K, K+1 ) / T
257: D = T*( AK*AKP1-ONE )
258: A( K, K ) = AKP1 / D
259: A( K+1, K+1 ) = AK / D
260: A( K, K+1 ) = -AKKP1 / D
261: *
262: * Compute columns K and K+1 of the inverse.
263: *
264: IF( K.GT.1 ) THEN
265: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
266: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
267: $ A( 1, K ), 1 )
268: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
269: $ K ), 1 ) )
270: A( K, K+1 ) = A( K, K+1 ) -
271: $ ZDOTC( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
272: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
273: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
274: $ A( 1, K+1 ), 1 )
275: A( K+1, K+1 ) = A( K+1, K+1 ) -
276: $ DBLE( ZDOTC( K-1, WORK, 1, A( 1, K+1 ),
277: $ 1 ) )
278: END IF
279: KSTEP = 2
280: END IF
281: *
282: IF( KSTEP.EQ.1 ) THEN
283: *
284: * Interchange rows and columns K and IPIV(K) in the leading
285: * submatrix A(1:k,1:k)
286: *
287: KP = IPIV( K )
288: IF( KP.NE.K ) THEN
289: *
290: IF( KP.GT.1 )
291: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
292: *
293: DO 40 J = KP + 1, K - 1
294: TEMP = DCONJG( A( J, K ) )
295: A( J, K ) = DCONJG( A( KP, J ) )
296: A( KP, J ) = TEMP
297: 40 CONTINUE
298: *
299: A( KP, K ) = DCONJG( A( KP, K ) )
300: *
301: TEMP = A( K, K )
302: A( K, K ) = A( KP, KP )
303: A( KP, KP ) = TEMP
304: END IF
305: ELSE
306: *
307: * Interchange rows and columns K and K+1 with -IPIV(K) and
308: * -IPIV(K+1) in the leading submatrix A(k+1:n,k+1:n)
309: *
310: * (1) Interchange rows and columns K and -IPIV(K)
311: *
312: KP = -IPIV( K )
313: IF( KP.NE.K ) THEN
314: *
315: IF( KP.GT.1 )
316: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
317: *
318: DO 50 J = KP + 1, K - 1
319: TEMP = DCONJG( A( J, K ) )
320: A( J, K ) = DCONJG( A( KP, J ) )
321: A( KP, J ) = TEMP
322: 50 CONTINUE
323: *
324: A( KP, K ) = DCONJG( A( KP, K ) )
325: *
326: TEMP = A( K, K )
327: A( K, K ) = A( KP, KP )
328: A( KP, KP ) = TEMP
329: *
330: TEMP = A( K, K+1 )
331: A( K, K+1 ) = A( KP, K+1 )
332: A( KP, K+1 ) = TEMP
333: END IF
334: *
335: * (2) Interchange rows and columns K+1 and -IPIV(K+1)
336: *
337: K = K + 1
338: KP = -IPIV( K )
339: IF( KP.NE.K ) THEN
340: *
341: IF( KP.GT.1 )
342: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
343: *
344: DO 60 J = KP + 1, K - 1
345: TEMP = DCONJG( A( J, K ) )
346: A( J, K ) = DCONJG( A( KP, J ) )
347: A( KP, J ) = TEMP
348: 60 CONTINUE
349: *
350: A( KP, K ) = DCONJG( A( KP, K ) )
351: *
352: TEMP = A( K, K )
353: A( K, K ) = A( KP, KP )
354: A( KP, KP ) = TEMP
355: END IF
356: END IF
357: *
358: K = K + 1
359: GO TO 30
360: 70 CONTINUE
361: *
362: ELSE
363: *
364: * Compute inv(A) from the factorization A = L*D*L**H.
365: *
366: * K is the main loop index, decreasing from N to 1 in steps of
367: * 1 or 2, depending on the size of the diagonal blocks.
368: *
369: K = N
370: 80 CONTINUE
371: *
372: * If K < 1, exit from loop.
373: *
374: IF( K.LT.1 )
375: $ GO TO 120
376: *
377: IF( IPIV( K ).GT.0 ) THEN
378: *
379: * 1 x 1 diagonal block
380: *
381: * Invert the diagonal block.
382: *
383: A( K, K ) = ONE / DBLE( A( K, K ) )
384: *
385: * Compute column K of the inverse.
386: *
387: IF( K.LT.N ) THEN
388: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
389: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
390: $ 1, CZERO, A( K+1, K ), 1 )
391: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
392: $ A( K+1, K ), 1 ) )
393: END IF
394: KSTEP = 1
395: ELSE
396: *
397: * 2 x 2 diagonal block
398: *
399: * Invert the diagonal block.
400: *
401: T = ABS( A( K, K-1 ) )
402: AK = DBLE( A( K-1, K-1 ) ) / T
403: AKP1 = DBLE( A( K, K ) ) / T
404: AKKP1 = A( K, K-1 ) / T
405: D = T*( AK*AKP1-ONE )
406: A( K-1, K-1 ) = AKP1 / D
407: A( K, K ) = AK / D
408: A( K, K-1 ) = -AKKP1 / D
409: *
410: * Compute columns K-1 and K of the inverse.
411: *
412: IF( K.LT.N ) THEN
413: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
414: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
415: $ 1, CZERO, A( K+1, K ), 1 )
416: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
417: $ A( K+1, K ), 1 ) )
418: A( K, K-1 ) = A( K, K-1 ) -
419: $ ZDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
420: $ 1 )
421: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
422: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
423: $ 1, CZERO, A( K+1, K-1 ), 1 )
424: A( K-1, K-1 ) = A( K-1, K-1 ) -
425: $ DBLE( ZDOTC( N-K, WORK, 1, A( K+1, K-1 ),
426: $ 1 ) )
427: END IF
428: KSTEP = 2
429: END IF
430: *
431: IF( KSTEP.EQ.1 ) THEN
432: *
433: * Interchange rows and columns K and IPIV(K) in the trailing
434: * submatrix A(k:n,k:n)
435: *
436: KP = IPIV( K )
437: IF( KP.NE.K ) THEN
438: *
439: IF( KP.LT.N )
440: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
441: *
442: DO 90 J = K + 1, KP - 1
443: TEMP = DCONJG( A( J, K ) )
444: A( J, K ) = DCONJG( A( KP, J ) )
445: A( KP, J ) = TEMP
446: 90 CONTINUE
447: *
448: A( KP, K ) = DCONJG( A( KP, K ) )
449: *
450: TEMP = A( K, K )
451: A( K, K ) = A( KP, KP )
452: A( KP, KP ) = TEMP
453: END IF
454: ELSE
455: *
456: * Interchange rows and columns K and K-1 with -IPIV(K) and
457: * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
458: *
459: * (1) Interchange rows and columns K and -IPIV(K)
460: *
461: KP = -IPIV( K )
462: IF( KP.NE.K ) THEN
463: *
464: IF( KP.LT.N )
465: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
466: *
467: DO 100 J = K + 1, KP - 1
468: TEMP = DCONJG( A( J, K ) )
469: A( J, K ) = DCONJG( A( KP, J ) )
470: A( KP, J ) = TEMP
471: 100 CONTINUE
472: *
473: A( KP, K ) = DCONJG( A( KP, K ) )
474: *
475: TEMP = A( K, K )
476: A( K, K ) = A( KP, KP )
477: A( KP, KP ) = TEMP
478: *
479: TEMP = A( K, K-1 )
480: A( K, K-1 ) = A( KP, K-1 )
481: A( KP, K-1 ) = TEMP
482: END IF
483: *
484: * (2) Interchange rows and columns K-1 and -IPIV(K-1)
485: *
486: K = K - 1
487: KP = -IPIV( K )
488: IF( KP.NE.K ) THEN
489: *
490: IF( KP.LT.N )
491: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
492: *
493: DO 110 J = K + 1, KP - 1
494: TEMP = DCONJG( A( J, K ) )
495: A( J, K ) = DCONJG( A( KP, J ) )
496: A( KP, J ) = TEMP
497: 110 CONTINUE
498: *
499: A( KP, K ) = DCONJG( A( KP, K ) )
500: *
501: TEMP = A( K, K )
502: A( K, K ) = A( KP, KP )
503: A( KP, KP ) = TEMP
504: END IF
505: END IF
506: *
507: K = K - 1
508: GO TO 80
509: 120 CONTINUE
510: END IF
511: *
512: RETURN
513: *
514: * End of ZHETRI_ROOK
515: *
516: 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