LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
csytri_rook.f
Go to the documentation of this file.
1*> \brief \b CSYTRI_ROOK
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytri_rook.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytri_rook.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytri_rook.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CSYTRI_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 A( LDA, * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CSYTRI_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 CSYTRF_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 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 CSYTRF_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 CSYTRF_ROOK.
86*> \endverbatim
87*>
88*> \param[out] WORK
89*> \verbatim
90*> WORK is COMPLEX 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 hetri_rook
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 csytri_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 A( LDA, * ), WORK( * )
141* ..
142*
143* =====================================================================
144*
145* .. Parameters ..
146 COMPLEX CONE, CZERO
147 parameter( cone = ( 1.0e+0, 0.0e+0 ),
148 \$ czero = ( 0.0e+0, 0.0e+0 ) )
149* ..
150* .. Local Scalars ..
151 LOGICAL UPPER
152 INTEGER K, KP, KSTEP
153 COMPLEX AK, AKKP1, AKP1, D, T, TEMP
154* ..
155* .. External Functions ..
156 LOGICAL LSAME
157 COMPLEX CDOTU
158 EXTERNAL lsame, cdotu
159* ..
160* .. External Subroutines ..
161 EXTERNAL ccopy, cswap, csymv, 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( 'CSYTRI_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 ccopy( k-1, a( 1, k ), 1, work, 1 )
237 CALL csymv( uplo, k-1, -cone, a, lda, work, 1, czero,
238 \$ a( 1, k ), 1 )
239 a( k, k ) = a( k, k ) - cdotu( 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 ccopy( k-1, a( 1, k ), 1, work, 1 )
262 CALL csymv( uplo, k-1, -cone, a, lda, work, 1, czero,
263 \$ a( 1, k ), 1 )
264 a( k, k ) = a( k, k ) - cdotu( k-1, work, 1, a( 1, k ),
265 \$ 1 )
266 a( k, k+1 ) = a( k, k+1 ) -
267 \$ cdotu( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
268 CALL ccopy( k-1, a( 1, k+1 ), 1, work, 1 )
269 CALL csymv( 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 \$ cdotu( 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 cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
286 CALL cswap( 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 cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
300 CALL cswap( 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 cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
315 CALL cswap( 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 ccopy( n-k, a( k+1, k ), 1, work, 1 )
353 CALL csymv( 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 ) - cdotu( 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 ccopy( n-k, a( k+1, k ), 1, work, 1 )
378 CALL csymv( 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 ) - cdotu( n-k, work, 1, a( k+1, k ),
381 \$ 1 )
382 a( k, k-1 ) = a( k, k-1 ) -
383 \$ cdotu( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
384 \$ 1 )
385 CALL ccopy( n-k, a( k+1, k-1 ), 1, work, 1 )
386 CALL csymv( 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 \$ cdotu( 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 cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
403 CALL cswap( 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 cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
417 CALL cswap( 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 cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
432 CALL cswap( 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 CSYTRI_ROOK
447*
448 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine csymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
CSYMV computes a matrix-vector product for a complex symmetric matrix.
Definition csymv.f:157
subroutine csytri_rook(uplo, n, a, lda, ipiv, work, info)
CSYTRI_ROOK
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81