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