1: *> \brief \b ZSYTRI_ROOK
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZSYTRI_ROOK + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytri_rook.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytri_rook.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytri_rook.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZSYTRI_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: *> ZSYTRI_ROOK computes the inverse of a complex symmetric
39: *> matrix A using the factorization A = U*D*U**T or A = L*D*L**T
40: *> computed by ZSYTRF_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**T;
52: *> = 'L': Lower triangular, form is A = L*D*L**T.
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 ZSYTRF_ROOK.
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 block structure of D
85: *> as determined by ZSYTRF_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 2011
111: *
112: *> \ingroup complex16SYcomputational
113: *
114: *> \par Contributors:
115: * ==================
116: *>
117: *> \verbatim
118: *>
119: *> November 2011, 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: *>
127: *> \endverbatim
128: *
129: * =====================================================================
130: SUBROUTINE ZSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
131: *
132: * -- LAPACK computational routine (version 3.4.0) --
133: * -- LAPACK is a software package provided by Univ. of Tennessee, --
134: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135: * November 2011
136: *
137: * .. Scalar Arguments ..
138: CHARACTER UPLO
139: INTEGER INFO, LDA, N
140: * ..
141: * .. Array Arguments ..
142: INTEGER IPIV( * )
143: COMPLEX*16 A( LDA, * ), WORK( * )
144: * ..
145: *
146: * =====================================================================
147: *
148: * .. Parameters ..
149: COMPLEX*16 CONE, CZERO
150: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
151: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
152: * ..
153: * .. Local Scalars ..
154: LOGICAL UPPER
155: INTEGER K, KP, KSTEP
156: COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP
157: * ..
158: * .. External Functions ..
159: LOGICAL LSAME
160: COMPLEX*16 ZDOTU
161: EXTERNAL LSAME, ZDOTU
162: * ..
163: * .. External Subroutines ..
164: EXTERNAL ZCOPY, ZSWAP, ZSYMV, XERBLA
165: * ..
166: * .. Intrinsic Functions ..
167: INTRINSIC MAX
168: * ..
169: * .. Executable Statements ..
170: *
171: * Test the input parameters.
172: *
173: INFO = 0
174: UPPER = LSAME( UPLO, 'U' )
175: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
176: INFO = -1
177: ELSE IF( N.LT.0 ) THEN
178: INFO = -2
179: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
180: INFO = -4
181: END IF
182: IF( INFO.NE.0 ) THEN
183: CALL XERBLA( 'ZSYTRI_ROOK', -INFO )
184: RETURN
185: END IF
186: *
187: * Quick return if possible
188: *
189: IF( N.EQ.0 )
190: $ RETURN
191: *
192: * Check that the diagonal matrix D is nonsingular.
193: *
194: IF( UPPER ) THEN
195: *
196: * Upper triangular storage: examine D from bottom to top
197: *
198: DO 10 INFO = N, 1, -1
199: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
200: $ RETURN
201: 10 CONTINUE
202: ELSE
203: *
204: * Lower triangular storage: examine D from top to bottom.
205: *
206: DO 20 INFO = 1, N
207: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
208: $ RETURN
209: 20 CONTINUE
210: END IF
211: INFO = 0
212: *
213: IF( UPPER ) THEN
214: *
215: * Compute inv(A) from the factorization A = U*D*U**T.
216: *
217: * K is the main loop index, increasing from 1 to N in steps of
218: * 1 or 2, depending on the size of the diagonal blocks.
219: *
220: K = 1
221: 30 CONTINUE
222: *
223: * If K > N, exit from loop.
224: *
225: IF( K.GT.N )
226: $ GO TO 40
227: *
228: IF( IPIV( K ).GT.0 ) THEN
229: *
230: * 1 x 1 diagonal block
231: *
232: * Invert the diagonal block.
233: *
234: A( K, K ) = CONE / A( K, K )
235: *
236: * Compute column K of the inverse.
237: *
238: IF( K.GT.1 ) THEN
239: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
240: CALL ZSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
241: $ A( 1, K ), 1 )
242: A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
243: $ 1 )
244: END IF
245: KSTEP = 1
246: ELSE
247: *
248: * 2 x 2 diagonal block
249: *
250: * Invert the diagonal block.
251: *
252: T = A( K, K+1 )
253: AK = A( K, K ) / T
254: AKP1 = A( K+1, K+1 ) / T
255: AKKP1 = A( K, K+1 ) / T
256: D = T*( AK*AKP1-CONE )
257: A( K, K ) = AKP1 / D
258: A( K+1, K+1 ) = AK / D
259: A( K, K+1 ) = -AKKP1 / D
260: *
261: * Compute columns K and K+1 of the inverse.
262: *
263: IF( K.GT.1 ) THEN
264: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
265: CALL ZSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
266: $ A( 1, K ), 1 )
267: A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
268: $ 1 )
269: A( K, K+1 ) = A( K, K+1 ) -
270: $ ZDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
271: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
272: CALL ZSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
273: $ A( 1, K+1 ), 1 )
274: A( K+1, K+1 ) = A( K+1, K+1 ) -
275: $ ZDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 )
276: END IF
277: KSTEP = 2
278: END IF
279: *
280: IF( KSTEP.EQ.1 ) THEN
281: *
282: * Interchange rows and columns K and IPIV(K) in the leading
283: * submatrix A(1:k+1,1:k+1)
284: *
285: KP = IPIV( K )
286: IF( KP.NE.K ) THEN
287: IF( KP.GT.1 )
288: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
289: CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
290: TEMP = A( K, K )
291: A( K, K ) = A( KP, KP )
292: A( KP, KP ) = TEMP
293: END IF
294: ELSE
295: *
296: * Interchange rows and columns K and K+1 with -IPIV(K) and
297: * -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
298: *
299: KP = -IPIV( K )
300: IF( KP.NE.K ) THEN
301: IF( KP.GT.1 )
302: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
303: CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
304: *
305: TEMP = A( K, K )
306: A( K, K ) = A( KP, KP )
307: A( KP, KP ) = TEMP
308: TEMP = A( K, K+1 )
309: A( K, K+1 ) = A( KP, K+1 )
310: A( KP, K+1 ) = TEMP
311: END IF
312: *
313: K = K + 1
314: KP = -IPIV( K )
315: IF( KP.NE.K ) THEN
316: IF( KP.GT.1 )
317: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
318: CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
319: TEMP = A( K, K )
320: A( K, K ) = A( KP, KP )
321: A( KP, KP ) = TEMP
322: END IF
323: END IF
324: *
325: K = K + 1
326: GO TO 30
327: 40 CONTINUE
328: *
329: ELSE
330: *
331: * Compute inv(A) from the factorization A = L*D*L**T.
332: *
333: * K is the main loop index, increasing from 1 to N in steps of
334: * 1 or 2, depending on the size of the diagonal blocks.
335: *
336: K = N
337: 50 CONTINUE
338: *
339: * If K < 1, exit from loop.
340: *
341: IF( K.LT.1 )
342: $ GO TO 60
343: *
344: IF( IPIV( K ).GT.0 ) THEN
345: *
346: * 1 x 1 diagonal block
347: *
348: * Invert the diagonal block.
349: *
350: A( K, K ) = CONE / A( K, K )
351: *
352: * Compute column K of the inverse.
353: *
354: IF( K.LT.N ) THEN
355: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
356: CALL ZSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
357: $ CZERO, A( K+1, K ), 1 )
358: A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
359: $ 1 )
360: END IF
361: KSTEP = 1
362: ELSE
363: *
364: * 2 x 2 diagonal block
365: *
366: * Invert the diagonal block.
367: *
368: T = A( K, K-1 )
369: AK = A( K-1, K-1 ) / T
370: AKP1 = A( K, K ) / T
371: AKKP1 = A( K, K-1 ) / T
372: D = T*( AK*AKP1-CONE )
373: A( K-1, K-1 ) = AKP1 / D
374: A( K, K ) = AK / D
375: A( K, K-1 ) = -AKKP1 / D
376: *
377: * Compute columns K-1 and K of the inverse.
378: *
379: IF( K.LT.N ) THEN
380: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
381: CALL ZSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
382: $ CZERO, A( K+1, K ), 1 )
383: A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
384: $ 1 )
385: A( K, K-1 ) = A( K, K-1 ) -
386: $ ZDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
387: $ 1 )
388: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
389: CALL ZSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
390: $ CZERO, A( K+1, K-1 ), 1 )
391: A( K-1, K-1 ) = A( K-1, K-1 ) -
392: $ ZDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 )
393: END IF
394: KSTEP = 2
395: END IF
396: *
397: IF( KSTEP.EQ.1 ) THEN
398: *
399: * Interchange rows and columns K and IPIV(K) in the trailing
400: * submatrix A(k-1:n,k-1:n)
401: *
402: KP = IPIV( K )
403: IF( KP.NE.K ) THEN
404: IF( KP.LT.N )
405: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
406: CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
407: TEMP = A( K, K )
408: A( K, K ) = A( KP, KP )
409: A( KP, KP ) = TEMP
410: END IF
411: ELSE
412: *
413: * Interchange rows and columns K and K-1 with -IPIV(K) and
414: * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
415: *
416: KP = -IPIV( K )
417: IF( KP.NE.K ) THEN
418: IF( KP.LT.N )
419: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
420: CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
421: *
422: TEMP = A( K, K )
423: A( K, K ) = A( KP, KP )
424: A( KP, KP ) = TEMP
425: TEMP = A( K, K-1 )
426: A( K, K-1 ) = A( KP, K-1 )
427: A( KP, K-1 ) = TEMP
428: END IF
429: *
430: K = K - 1
431: KP = -IPIV( K )
432: IF( KP.NE.K ) THEN
433: IF( KP.LT.N )
434: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
435: CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
436: TEMP = A( K, K )
437: A( K, K ) = A( KP, KP )
438: A( KP, KP ) = TEMP
439: END IF
440: END IF
441: *
442: K = K - 1
443: GO TO 50
444: 60 CONTINUE
445: END IF
446: *
447: RETURN
448: *
449: * End of ZSYTRI_ROOK
450: *
451: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsytri_rook.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack