LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
csytri.f
Go to the documentation of this file.
1*> \brief \b CSYTRI
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.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CSYTRI( 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 computes the inverse of a complex symmetric indefinite matrix
39*> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
40*> CSYTRF.
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.
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.
86*> \endverbatim
87*>
88*> \param[out] WORK
89*> \verbatim
90*> WORK is COMPLEX array, dimension (2*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
111*
112* =====================================================================
113 SUBROUTINE csytri( UPLO, N, A, LDA, IPIV, WORK, INFO )
114*
115* -- LAPACK computational routine --
116* -- LAPACK is a software package provided by Univ. of Tennessee, --
117* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
118*
119* .. Scalar Arguments ..
120 CHARACTER UPLO
121 INTEGER INFO, LDA, N
122* ..
123* .. Array Arguments ..
124 INTEGER IPIV( * )
125 COMPLEX A( LDA, * ), WORK( * )
126* ..
127*
128* =====================================================================
129*
130* .. Parameters ..
131 COMPLEX ONE, ZERO
132 parameter( one = ( 1.0e+0, 0.0e+0 ),
133 \$ zero = ( 0.0e+0, 0.0e+0 ) )
134* ..
135* .. Local Scalars ..
136 LOGICAL UPPER
137 INTEGER K, KP, KSTEP
138 COMPLEX AK, AKKP1, AKP1, D, T, TEMP
139* ..
140* .. External Functions ..
141 LOGICAL LSAME
142 COMPLEX CDOTU
143 EXTERNAL lsame, cdotu
144* ..
145* .. External Subroutines ..
146 EXTERNAL ccopy, cswap, csymv, xerbla
147* ..
148* .. Intrinsic Functions ..
149 INTRINSIC abs, max
150* ..
151* .. Executable Statements ..
152*
153* Test the input parameters.
154*
155 info = 0
156 upper = lsame( uplo, 'U' )
157 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
158 info = -1
159 ELSE IF( n.LT.0 ) THEN
160 info = -2
161 ELSE IF( lda.LT.max( 1, n ) ) THEN
162 info = -4
163 END IF
164 IF( info.NE.0 ) THEN
165 CALL xerbla( 'CSYTRI', -info )
166 RETURN
167 END IF
168*
169* Quick return if possible
170*
171 IF( n.EQ.0 )
172 \$ RETURN
173*
174* Check that the diagonal matrix D is nonsingular.
175*
176 IF( upper ) THEN
177*
178* Upper triangular storage: examine D from bottom to top
179*
180 DO 10 info = n, 1, -1
181 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
182 \$ RETURN
183 10 CONTINUE
184 ELSE
185*
186* Lower triangular storage: examine D from top to bottom.
187*
188 DO 20 info = 1, n
189 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
190 \$ RETURN
191 20 CONTINUE
192 END IF
193 info = 0
194*
195 IF( upper ) THEN
196*
197* Compute inv(A) from the factorization A = U*D*U**T.
198*
199* K is the main loop index, increasing from 1 to N in steps of
200* 1 or 2, depending on the size of the diagonal blocks.
201*
202 k = 1
203 30 CONTINUE
204*
205* If K > N, exit from loop.
206*
207 IF( k.GT.n )
208 \$ GO TO 40
209*
210 IF( ipiv( k ).GT.0 ) THEN
211*
212* 1 x 1 diagonal block
213*
214* Invert the diagonal block.
215*
216 a( k, k ) = one / a( k, k )
217*
218* Compute column K of the inverse.
219*
220 IF( k.GT.1 ) THEN
221 CALL ccopy( k-1, a( 1, k ), 1, work, 1 )
222 CALL csymv( uplo, k-1, -one, a, lda, work, 1, zero,
223 \$ a( 1, k ), 1 )
224 a( k, k ) = a( k, k ) - cdotu( k-1, work, 1, a( 1, k ),
225 \$ 1 )
226 END IF
227 kstep = 1
228 ELSE
229*
230* 2 x 2 diagonal block
231*
232* Invert the diagonal block.
233*
234 t = a( k, k+1 )
235 ak = a( k, k ) / t
236 akp1 = a( k+1, k+1 ) / t
237 akkp1 = a( k, k+1 ) / t
238 d = t*( ak*akp1-one )
239 a( k, k ) = akp1 / d
240 a( k+1, k+1 ) = ak / d
241 a( k, k+1 ) = -akkp1 / d
242*
243* Compute columns K and K+1 of the inverse.
244*
245 IF( k.GT.1 ) THEN
246 CALL ccopy( k-1, a( 1, k ), 1, work, 1 )
247 CALL csymv( uplo, k-1, -one, a, lda, work, 1, zero,
248 \$ a( 1, k ), 1 )
249 a( k, k ) = a( k, k ) - cdotu( k-1, work, 1, a( 1, k ),
250 \$ 1 )
251 a( k, k+1 ) = a( k, k+1 ) -
252 \$ cdotu( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
253 CALL ccopy( k-1, a( 1, k+1 ), 1, work, 1 )
254 CALL csymv( uplo, k-1, -one, a, lda, work, 1, zero,
255 \$ a( 1, k+1 ), 1 )
256 a( k+1, k+1 ) = a( k+1, k+1 ) -
257 \$ cdotu( k-1, work, 1, a( 1, k+1 ), 1 )
258 END IF
259 kstep = 2
260 END IF
261*
262 kp = abs( ipiv( k ) )
263 IF( kp.NE.k ) THEN
264*
265* Interchange rows and columns K and KP in the leading
266* submatrix A(1:k+1,1:k+1)
267*
268 CALL cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
269 CALL cswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
270 temp = a( k, k )
271 a( k, k ) = a( kp, kp )
272 a( kp, kp ) = temp
273 IF( kstep.EQ.2 ) THEN
274 temp = a( k, k+1 )
275 a( k, k+1 ) = a( kp, k+1 )
276 a( kp, k+1 ) = temp
277 END IF
278 END IF
279*
280 k = k + kstep
281 GO TO 30
282 40 CONTINUE
283*
284 ELSE
285*
286* Compute inv(A) from the factorization A = L*D*L**T.
287*
288* K is the main loop index, increasing from 1 to N in steps of
289* 1 or 2, depending on the size of the diagonal blocks.
290*
291 k = n
292 50 CONTINUE
293*
294* If K < 1, exit from loop.
295*
296 IF( k.LT.1 )
297 \$ GO TO 60
298*
299 IF( ipiv( k ).GT.0 ) THEN
300*
301* 1 x 1 diagonal block
302*
303* Invert the diagonal block.
304*
305 a( k, k ) = one / a( k, k )
306*
307* Compute column K of the inverse.
308*
309 IF( k.LT.n ) THEN
310 CALL ccopy( n-k, a( k+1, k ), 1, work, 1 )
311 CALL csymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
312 \$ zero, a( k+1, k ), 1 )
313 a( k, k ) = a( k, k ) - cdotu( n-k, work, 1, a( k+1, k ),
314 \$ 1 )
315 END IF
316 kstep = 1
317 ELSE
318*
319* 2 x 2 diagonal block
320*
321* Invert the diagonal block.
322*
323 t = a( k, k-1 )
324 ak = a( k-1, k-1 ) / t
325 akp1 = a( k, k ) / t
326 akkp1 = a( k, k-1 ) / t
327 d = t*( ak*akp1-one )
328 a( k-1, k-1 ) = akp1 / d
329 a( k, k ) = ak / d
330 a( k, k-1 ) = -akkp1 / d
331*
332* Compute columns K-1 and K of the inverse.
333*
334 IF( k.LT.n ) THEN
335 CALL ccopy( n-k, a( k+1, k ), 1, work, 1 )
336 CALL csymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
337 \$ zero, a( k+1, k ), 1 )
338 a( k, k ) = a( k, k ) - cdotu( n-k, work, 1, a( k+1, k ),
339 \$ 1 )
340 a( k, k-1 ) = a( k, k-1 ) -
341 \$ cdotu( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
342 \$ 1 )
343 CALL ccopy( n-k, a( k+1, k-1 ), 1, work, 1 )
344 CALL csymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
345 \$ zero, a( k+1, k-1 ), 1 )
346 a( k-1, k-1 ) = a( k-1, k-1 ) -
347 \$ cdotu( n-k, work, 1, a( k+1, k-1 ), 1 )
348 END IF
349 kstep = 2
350 END IF
351*
352 kp = abs( ipiv( k ) )
353 IF( kp.NE.k ) THEN
354*
355* Interchange rows and columns K and KP in the trailing
356* submatrix A(k-1:n,k-1:n)
357*
358 IF( kp.LT.n )
359 \$ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
360 CALL cswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
361 temp = a( k, k )
362 a( k, k ) = a( kp, kp )
363 a( kp, kp ) = temp
364 IF( kstep.EQ.2 ) THEN
365 temp = a( k, k-1 )
366 a( k, k-1 ) = a( kp, k-1 )
367 a( kp, k-1 ) = temp
368 END IF
369 END IF
370*
371 k = k - kstep
372 GO TO 50
373 60 CONTINUE
374 END IF
375*
376 RETURN
377*
378* End of CSYTRI
379*
380 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(uplo, n, a, lda, ipiv, work, info)
CSYTRI
Definition csytri.f:114
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81